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and Kelley [ 3 ] . However, these papers were not concerned with 

promoting the basic understanding of the traction between elastohydro- 
dynamic contacts. 

Crook [ 4 J used two kinds of rolling disk machines in measuring 

the friction in a line contact as a function of sliding speed. In 
the region of small sliding speeds, he used the four-disk machine, a 
center disk surrounded by three equally spaced outer disks, shown in 
Figure 1.1. The center disk is free-floating and the measured torque 
does not contain any extraneous torque from the supporting bearings. 

For this reason, the four-disk machine gives very accurate frictional 
torque measurements at small sliding speeds. The four-disk machine 
is not suitable in the region of high slips, however, since it cannot 
maintain a stable sliding speed. For high sliding speeds, Crook used 
the two-disk machine shown in Figure 1.2, where the rotations of both 
disks are controlled by variable speed motors. Thus, Crook was able 
to measure the friction characteristics throughout the entire range 
of sliding speeds, using the four -disk machine in the low slip region 
and the two-disk machine in the high sliding speed region. 

Crook found a profound influence of rolling speed upon the 
frictional torque in the low slip region. In this region, the slope 
of the traction versus slip curve is equal to the "effective viscosity" 
divided by the oil film thickness. Therefore, the effective viscosity 
may be evaluated by measuring the slope of the traction curve and 
calculating the oil film thickness from existing elastohydrodynamic 
theory. If the thermal effects and the non-Newtonian effects of the 
lubricant were both absent in this region, the effective viscosity 
would not be a function of rolling speed. However, this condition was 
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Figure 1.2. Crook’s two-disk machine. A and B disks; C and D 

swinging arms; E axle; F and G loading cables; H spring 
beam; I dial gauge. Figure from Crook [ 4 ] . 
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not found in Crook's experimental results. On the contrary, he found 
a marked influence of the rolling speed on the effective viscosity of 
the lubricant that does not appear to be due to thermal effects only. 
Crook speculated that it was the viscoelastic effect of the lubricant 
which prevented it from reaching the static viscosity in the short 
time interval as it passes through the contact zone. 

Crook was able to extend the friction data in the high slip 

region, with his two-disk machine, for loads ranging from 7.5 to 20 x 10^ 
2 

dynes/cm and rolling speeds from 400 cm/sec to 1200 cm/sec. All the 
friction curves show the same basic trend which is characterized by 
an ascending portion at ‘small sliding speeds and a descending friction 
at high sliding speeds. An increase in load does not change the basic 
characteristics of the friction curve, but does increase the level of the 
friction force. Similarly, Crook found that an increase in the rolling 
speed decreases the friction level. 

Crook also attempted to predict the friction analytically by a 
simplified thermal friction theory based on the following four as- 
sumptions: the film thickness within the contact zone is uniform; 

the pressure distribution in the contact region is Hertzian; the heat 
carried away by the lubricant due to convection may be neglected; and 
the temperature rise on the surface of the disk may also be neglected. 
Using this simplified theory for a Newtonian lubricant. Crook was able 
to calculate the coefficient of friction or the effective viscosity 
as a function of sliding s.peed. However, he could not predict the 
sharp reduction of the effective viscosity at small sliding speeds. 

He concluded that the friction force at small sliding speeds cannot 
be accurately predicted by considering the thermal effects only. 
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Cheng [ 5 ] employed his full e lastohydrodynamic theory in 

calculating the friction for the conditions corresponding to those 
used in Crook's experiments. The temperature calculations are based 
on the finite difference solution of the energy equation and are free 
from all the assumptions made earlier by Crook. It is seen in Figure 
1.3 that even with this refined thermal analysis there still exists a 
large discrepancy in the low slip region. This strengthens Crook's 
argument that the thermal effects alone cannot account for the sharp 
reduction of effective viscosity in the low slip region. 

Bellj'Kannel and Allen [ 6 ] developed an approximate analysis 

to predict the temperature rise in the lubricant film at low sliding 

speeds. Their analysis included the heat due to convection and the 

heat generation due to the compression of the lubricant. They also 

concluded that the temperature effects are too small to account for 

the loss of effective viscosity at low sliding speeds. In addition 

to the thermal theory, they developed a non-Newtonian friction theory 

using a rheological model proposed by Ree and Eyring [ 7 ] . The 

results of this analysis indicate that drastic reductions of friction 

can exist if the lubricant viscosity is shear rate-dependent according 

to Ree-Eyring. However, in all their calculated data, the friction 

force was found to be dependent upon 1/h as the rolling speed is varied, 

whereas all the experimental data gathered thus far has shown the 

proportionality to be far greater than 1/h and in most cases more nearly 

2 

proportional to 1/h . Thus, the inclusion of the Ree-Eyring model 
alone in the friction analysis would not be able to predict a suf- 
ficient reduction of friction at low rolling speeds. 

Smith [ 8 J employed the rolling contact machine shown in Figure 1.4 
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Figure 1.3. Comparison of Cheng's therma 1 . theory with Crook's experiments. Curve from 
McGrew, Gu, Cheng and Murray [ 9 ] . 




Figure 1.4. Smith’s disk machine. 

a Cylindrical roller d Pivots g Strain gauge dynamometer 

b Spherical roller e Motor 
c Bearings f Normal load 

Figure from Smith [ 10 ] * 
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to measure the friction between two rollers whose axes are skewed at 
an arbitrary angle. With this skewed arrangement, he was able to 
measure the friction force due to the sliding velocity component. The 
resulting friction versus sliding speed curves show trends similar to 
those observed by Crook. Smith divided these curves into several 
regions. He believed a Newtonian isothermal friction theory is ap- 
plicable in the region where the friction varies proportionally with 
the sliding speed. In the ascending portion of the friction curves, 
where the friction force increases with sliding speed in a non-linear 
fashion, Smith believed that the non-linearity is due to the non- 
Newtonian behavior of the lubricant. He postulated that there is also 
a region in which a shear plane exists at the center of the lubricant 
film, such that the lubricant behaves like two solid layers sliding 
over each other at the shear plane. He further stated that the re- 
sistance to sliding at the shear plane is dependent upon the shear 
plane temperature and the hydrostatic stress in the lubricant. Finally, 
he defined a region where seizure would take place. 

A more comprehensive study of friction covering a wide range of 
loads, rolling speeds and sliding speeds was carried out more recently 
by Johnson and Cameron [ 11 ] with a two-disk machine. The maximum 

Hertzian pressure was varied from 62,000 psi to 243,000 psi; the rolling 
speed was varied from 8 in/sec to 260 in/sec; and the oil inlet tem- 
perature was varied from 30 °C to 90 °C. The most striking feature 
of Johnson and Cameron's data is that there exists a ceiling to all 
the experimental traction coefficients which cannot be exceeded no 
matter how the load and the rolling speed are varied. They also took 
extensive data in the low slip region, and from the slope of the 
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traction versus slip curve were able to calculate the effective vis- 
cosity as a function of rolling speed. Johnson and Cameron furnished 
more convincing evidence that the thermal effect s.: a lone cannot account 
for the experimentally measured friction, and that a Smith-type limiting 
shear stress is dependent only on the shear plane temperature and 
pressure . 

Jeffris and Johnson [12 ] investigated the effect of surface 

roughness upon friction between two lubricated rollers. They concluded 
that the measured coefficient of friction showed remarkably little 
variation throughout the whole range of experimental conditions for 
Hertzian pressures in excess of 175 kpsi. At lower Hertzian pressures, 
the surface roughness effect can be substantial. 

A rather interesting qualitative explanation of the velocity, 
rate of shear, viscosity and temperature variations across the film 
thickness of an elastohydrodynamic contact was offered by Plint [ 13 J 

He postulated that at the entrance of the contact zone, the rate of 
shear, viscosity and temperature are constant across the film thick- 
ness and the velocity profile is linear. As the thermal effects take 
over, the temperature at the mid-film increases and the viscosity is 
at a minimum at this position. This thermal effect causes the velocity 
profile to distort into a cusp such that the rate of shear becomes a 
maximum at the mid-film. A ceiling curve similar to that of Johnson 
and Cameron's was also found in Plint's experimental friction data. 

Dowson and Holmes [ 14 ] modified Crook's four-disk machine 

and investigated the effect of surface quality upon the traction char- 
acteristics of rolling and sliding contacts. They showed that the 
friction initially decreases with surface roughness, reaches a minimum, 
and then increases steadily with surface roughness. Unlike Jefferis 
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and Johnson’s conclusion on the effect of surface roughness on friction, 
Dowson and Holmes found that the influence of surface quality is quite 
pronounced. However, these two results may not be in direct contra- 
diction since the loads used by Dowson and Holmes were much smaller 
than those used in Jefferis and Johnson’s experiments. 

Recently, Dyson [ 15 j has made a pioneering study analyzing 
the frictional force in an elastohydrodynamic contact by considering 

, t 

the lubricant as a viscoelastic liquid. He simplified his analysis 
by dividing the friction versus sliding speed curve into three regions, 
as shown in Figure 1.5: the linear region, where the frictional force 

varies linearly with the sliding speed; the non-linear ascending region, 
where the slope of the friction curve reduces rapidly as the sliding 
speed increases; and the thermal region, where the frictional force 
decreases with the sliding speed. The results of this study are most 
encouraging and have inspired the author's investigation of the rhe- 
ological effects in an elastohydrodynamic lubricated contact. 

The friction analysis presented in this thesis describes the 
rheological behavior of the lubricant in an elastohydrodynamic concen- 
trated contact in terms of two viscoelastic models. These models 
represent the separate effects of non-Newtonian behavior and the 
transient response of the fluid. 

A unified description of the non-Newtonian shear rate dependence 
of the viscosity is presented in Chapter II as a new hyperbolic liquid 
model. The hyperbolic model is based upon a shear viscoelastic liquid 
model with the addition of a limiting value of shear stress. The limit 
ing shear stress is related to the high frequency limiting shear modulus 
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frictional traction 



Figure 1.5. Friction versus sliding speed curve. Curve from 
Dyson [15 ] . 
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of the lubricant , as proposed by Dyson [ 15 ] . 

The transient response of the viscosity, following the rapid 
pressure rise encountered in the contact, is described in Chapter III 
by a compressions 1 viscoelastic model of the volume response of a 
liquid to an applied pressure step. Kovacs [ 16 J first investigated 
this non-linear model for the volume creep of polymer melts. 

■ The governing equations, the fluid property functions and the 
technique used to calculate the tractive force transmitted during 
sliding between the two surfaces of a rolling disk machine are developed 
in Chapter IV. The experimental investigation is detailed in Chapter V 
and the analytical and experimental results are discussed and corre- 
lated in Chapter VI. 



CHAPTER II 


NON-LINEAR SHEAR STRESS -STRAIN RELATION 

A friction analysis based upon a Newtonian lubricant having a 
viscosity varying with the statically measured pressure and tempera- 
ture can yield a frictional coefficient far greater than those measured. 
There is little doubt that the fluid ceases to be Newtonian. Thus, a 
realistic friction analysis must consider a non-linear relationship 
between shear stress and shear rate. 

A major difficulty in predicting the friction for elastohydro- 
dynamic lubrication is the lack of data available for the physical 
properties of the lubricant at the extreme values of pressure, tempera- 
ture and shear rate encountered in the concentrated contact. It is 
very difficult to make direct measurements of shear stress versus 
shear rate in continuous shear under EHD conditions. Therefore, the 
physical property data must come from other fields. One source is 
the study of supercooled liquids under oscillatory shear. A corre- 
lation between this data and the behavior under conditions of continuous 
shear, as well as the restrictions of such a correlation, are discussed 
in this chapter. A hyperbolic shear stress-strain function is then 
proposed as a useful non-linear model for use under elastohydrodynamic 
conditions . 

2.1 Shear Viscoelasticity 

The phenomenological theory of viscoelasticity attempts to 
describe the mechanical behavior of a material in terms of time-de- 
pendent, or frequency-dependent, functions which relate the stress in 
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the material to the deformation. Classical elasticity theory for 
solids is based on Hooke's Law which requires the stress to be di- 
rectly proportional to the instantaneous strain but independent of the 
rate of strain. Classical hydrodynamic theory is based on Newton's 
Law. Newton’s Law states that the steady-state shear stress in a 
liquid is directly proportional to the instantaneous rate of strain 
but independent of the strain itself. Many materials closely follow 
the behavior specified by these laws. It is often impossible, however, 
to characterize a material by either of the two classical types of 
behavior. Substances which exhibit both solid-like and liquid-like 
properties show viscoelastic behavior. 

The term viscoelastic is used to describe the properties of any 
material which is able to store energy in elastic deformation and 
dissipate energy as heat. If the strain and strain rate are kept 
sufficient small, so that in a given experiment the ratio of stress 
to strain is a function of time only and independent of the stress 
level, the material shows linear viscoelastic behavior. Most of the 
physical properties of viscoelastic materials have been determined 
by oscillatory shear experiments. Linear viscoelastic behavior is 
easily obtained in these experiments since the amplitude of deform- 
ation is extremely small. 

The work of Gross [ 17 ] and Alfrey [ 18 J are examples of the 
large literature concerning the mathematical aspects of the phenom- 
enological theory of linear viscoelasticity. It is more appropriate 
here to develop the subject in terms of simple mechanical models. 

The model approach is easier to understand and more c losely vrelated 
to the physical behavior of the materials. 
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2.2 Viscoelastic Functions 


In most of the high frequency techniques used for measuring the 
viscoelastic properties of liquids, a plane shear wave is propagated 
through the liquid. The shear stress t and the shear strain y are 
related by a complex quantity, the shear modulus 

G* = — (2.1) 

Y 

In a Hookean solid, the shear modulus is a real quantity since 
the stress varies in phase with the strain. In a Newtonian liquid, 
the stress is 90° out of phase with the strain. In the latter case, 
the shear modulus is an imaginary quantity and is determined from 
Newton's Law. The strain rate is represented by 

Y = “ iu>Y (2.2) 

and therefore the stress is calculated as 

T = r|Y = imyri (2.3) 

The shear modulus is now calculated by its definition, equation (2.1). 

G = iurn (2.4) 

For a viscoelastic' material, the stress and strain differ by a phase 
angle between 0° and 90°. Therefore, the frequency-dependent shear 
modulus is a complex quantity with both real and imaginary components, 
as represented by 

G ( ict>) = G'(ou) + i G"(cu) (2.5) 

The shear modulus will not have the simple form given in equa- 
tion (2.4) for a Newtonian liquid except at low frequencies where 
sufficient time is available during each stress cycle for viscous flow 
to occur. At higher frequencies, the time required for molecular 
translation becomes comparable with the period of the stress cycle and 
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the liquid exhibits a shear rigidity. At sufficiently high frequencies, 
the behavior will be purely elastic. There is no molecular transition 
during each cycle and, consequently, the energy loss due to viscous 
flow is negligible. Under these conditions, the liquid behaves like 
a glass. 

The real component of the complex modulus G*, the ratio of the 
stress in phase with the strain to the strain, is called the storage 
modulus because of its association with the storage and release of 
elastic energy. The imaginary component G", the ratio of the stress 
90° out of phase with the strain to the strain, is called the loss 
modulus because of its association with the dissipation of energy as 
heat by viscous flow. 

The modulus components for a liquid have the following limits. 

At low frequencies where the behavior is purely viscous, or Newtonian: 
Lim G'(cu) = 0 (2.6) 

(jj -* 0 

Lim G"(ou) — u)T] (2.7) 

u) 0 

At high frequencies where the behavior is purely elastic: 

Lim G'(m) = G (2.8) 

CO 

U) -* » 

Lim G"(uj) = 0 (2.9) 

0) -* eo 

where G is the limiting elastic modulus. 

03 

☆ 

The shear mechanical impedance Z , defined as the ratio of shear 
stress to particle velocity, is the quantity most easily measured in 
the oscillatory experiments. It is mathematically related to the 
shear modulus by the equations governing shear wave propagation through 
a liquid medium. Barlow and Lamb [ 19 ] show this relationship to be 
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( 2 . 10 ) 


"k 0 k 

(Z ) = pG (la>) 

where p is the density of the liquid. 

k 

For a Newtonian liquid, where G is given by equation (2.4), 
the real and imaginary components of the shear mechanical impedance 
are given by 


z* = Z' + i Z" = (1 + i) J ^ 


( 2 . 11 ) 


Equation (2.10) allows the components of the shear modulus to 
be calculated from the experimentally measured components of the shear 
mechanical impedance as follows: 


G'(uj) 

G"(oo) 


(Z») 2 - (Z") 2 
P 

2 Z' Z" 

P 


( 2 . 12 ) 

(2.13) 


The liquid properties may be alternatively represented by a 
complex viscosity defined by 

* 1 

ri (iuu) = ri ’ (a>) - i-n"(a)) = ~ (2.14) 
The definition requires 

, (u))= £lM (2.15) 

co 

and 

*"(«) = (2.16) 

00 

The low frequency limit of the dynamic viscosity r)' is r) , the steady 
flow Newtonian viscosity. 

k 

The complex compliance of J is the inverse of the complex 
modulus . 

J*(iu>) = J*(oo) - i J"(uo) = 1 (2.17) 

G (iuj) 
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It follows that 


and 


J« = 


J" = 


G« 


[(G») 2 + (G") 2 ] 


G" 


(2.18) 


(2.19) 


[(G*) 2 + (G") 2 ] 

The real component J’ is the storage compliance and J" is the loss 
compliance . 

* 

The inverse of the complex viscosity is the complex fluidity p, . 

( 2 . 20 ) 


p (iu>) = p 1 (a>) + ip"(u>) = ■ 


ri (itn) 

The definitions and interrelations of the viscoelastic functions 
are summarized in Table 2.1. 


2,3 The Maxwell Model 

It is often convenient to visualize the behavior of a complex 
material in terms of models. The basic mechanical model elements are 
a coiled spring to represent Hookean elastic deformation and a dashpot 
to represent Newtonian viscous flow. Extension of the elements is 
analogous to shear strain and the associated force is analogous to the 
shear stress. 

The combination of a spring in series with a dashpot was studied 
by Maxwell [ 20 ] . This simple model, shown in Figure 2.1, exhibits 
both viscous and elastic behavior. Viscous flow in the dashpot with 
negligible extension of the spring takes place if the extension rate 
is small. If the model is rapidly extended and immediately released, 
the deformation is purely elastic since sufficient time is not avail- 
able for flow to occur in the' dashpot. Between these extremes, the 
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Table 2.1 ' 


VISCOELASTIC FUNCTIONS IN OSCILLATORY SHEAR 


T = shear stress acting in x-direction on x-z plane 
(• = particle displacement in x-direction 

Y = = shear strain 

• 

? = dE/^t = particle velocity in x-direction 

Y “ 3Y/d c = d|/dy = shear rate, rate of strain 
0 ) = angular frequency 


DEFINITIONS : 


Complex shear modulus: 

Complex mechanical impedance: 
Complex viscosity: 

Complex shear compliance: 
Complex fluidity: 


G (iu)) = t/y = G '(u)) + i 
2 (in>) = -t /? = z ’(u)) + 
T) (ito) = t/y = n'Ou) - i 

J*(iu>) - Y / t = - i 

p, (io>) = Y /t = vjl * (eo) + i 


INTERRELATIONS : 


* 

G 


(irn) 


— Z — — = iu)ri(Ljo) = 

■ P 


J (i«)) 


10 ) 

i(i(J)) 


G"(co) 

- Z"(co) 

r|"(o)) 

J"(u>) 

p"(u>) 
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behavior will be a combination of both the elastic and the viscous modes. 


~VSlQ£2jL2JiJlJL/~ 


Figure 2.1. The Maxwell element. The spring corresponds to a shear 
modulus G and the dashpot corresponds to a viscosity n . 


The basic equations of motion for the components of the model are: 

T = T)Y n (2.21) 

for the dashpot; and 


T = Gy. 


H 


( 2 . 22 ) 


for the spring, where 

T = the applied stress 
y N = the rate of extension of the dashpot 

y^ = the extension of the spring 


The rate of extension of the spring is = t/G and the total rate of 
extension is then 


• • i • T . 

= Y N + V H = j- + 


(2.23) 


or 


T + Q T - TlY 


(2.24) 


Equation (2.24) is the constitutive equation of the Maxwell element. 

The ratio r|/G has the dimensions of time and is called the Maxwell 
relation time \ . 

X = q (2.25) 
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For sinusoidal variations of stress and strain of frequency oo 


equation (2.24) becomes 

t + iuA.T = iojriY (2.26) 

The complex shear modulus is evaluated from equation (2.1) as 


io£L 


G (*»> = TTi 




Rationalizing this expression yields 
G*(icu) - 

1 + a) X 


(2.27) 


(2.28) 


and substituting for r) from equation (2.25) gives the final form of 
the complex shear modulus 

G*(iou) = G • ^ (2.29) 
1 + 0) X 


The storage modulus is 

2,2 

G*(cu) = G • (2.30) 

1 + «' V 

which reduces to G’((u) = G in the limit as m » ; but this limiting 
value has been defined as G . Thus, the spring in the Maxwell element 

CO 

corresponds to the instantaneous or limiting high frequency shear 
modulus of a liquid. The loss modulus is given by 

G"(u>) = G , (2.31) 

” 1 + u) x 

which in the limit as cu 0 becomes G"(u)) “ G • ooX = con • The dashpot 
of the Maxwell element therefore corresponds to the steady flow vis- 
cosity of a liquid. In normalized form, the variation with frequency 
of the modulus components and the dynamic viscosity is given by equa- 
tions (2.32), (2.33) and (2.34). 



(2.33) 


G"((b) _ mX 


G 

CO 


1 + 


2 2 
w \ 


El _ G " fa) 

T) CUT] 


, ^ 2 2 
1 + u> X 


(2.34) 


The complex compliance of J for the Maxwell element is given by the 
simple expression ' 


* 

J (in)) 


X _ 1 + iuA _ 

T iiUT) G 

00 



(2.35) 


The frequency variation of the modulus components and the dynamic 
viscosity for the Maxweil element are shown in Figure 2.2. 

Gruber and Litovitz [ 21 ] have postulated that a Maxwell element 
can predict the behavior of certain liquids. The viscosity of these 
liquids is governed primarily by the energy required for a molecule to 
surmount the potential barrier due to interaction with its nearest 
neighbors, and jump from one site in the liquid to another. The steady 
flow viscosity of such a liquid is given by the Arrhenius equation: 

In n = A + B/T (2.36) 

where T is the absolute temperature. 

For liquids which have viscosities above about 0.1 poise, however, 
the viscosity is primarily a function of the relative availability of 
free volume as described by Barlow, Lamb and Matheson [ 22 ] . There- 
fore, the Arrhenius viscosity-temperature relation and the Maxwell 
description of viscoelastic relaxation are not adequate governing equa- 
tions. The lubricant in an EHD concentrated contact is in a state 
where the viscosity is limited by available free volume, and, therefore, 
another viscoelastic liquid model proposed for such liquids by Barlow, 
Erginsav and Lamb [ 23 ] must be investigated. 
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MAXWELL MODEL: 
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Figure 2.2. The frequency variation 
the Maxwell element. 



2.4 The B. E. L. Liquid Model 


Barlow, Lamb and Matheson f 22 ]have shown that liquids having 
viscosities above about 0.1 poise obey the Doolittle free-volume T 24 j 
equation: 

v 

In t) = A + B — (2.37) 

V f 

where T) = viscosity 

v q = occupied volume 

v^ = free volume 

A,B = constants of a given liquid 

The specific volume v = v q + v^ and the density is a linear function 

of temperature. Therefore, equation (2.37) becomes 

In n = A* + B'/(T-T ) (2.38) 

' o 

where r| = viscosity at temperature T 

T = absolute temperature 

T q = reference temperature, at which there would be 
no free volume 

A*,B' = constants for a given liquid 
Barlow, Lamb, Matheson, Padmini and Richter f 25 ] and Barlow, 
Erginsav and Lamb T 23 ] have demonstrated that the viscoelastic 
properties for a large number and wide variety of liquids, which obey 

the Doolittle viscosity-temperature relation, can be represented by 

JL i 

two standard curves: Z'/(qG^) 2 and Z"/(pG^) 2 versus log^gfori/G^) . 

Figure 2.3 shows that the experimental results for many liquids are 
indistinguishable when plotted in this manner. This suggests a simple 
underlying phenomenological explanation. Barlow, Erginsav and Lamb 
propose a new liquid model consisting of the parallel combination of 
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the shear mechanical impedances for a Newtonian liquid and a Hookean 
solid. The shear mechanical impedances result from equations (2.8), 
(2.10) and (2.11). Thus, 

i 

Z N =(l+i )(^ 2 ( 2 . 39 ) 

for a Newtonian liquid and 

Z R = (pG (2.40) 


for a Hookean solid. 

Accordingly, the components of the shear mechanical impedance are 
given by: 


(pG <a ) 2 ((nn/ 2G oo )2 [ 

1 + (WGJ 2 j 


1 + (C0T1/2GJ2 

2 + 


(2.41) 


Z" = 


(pG {o )2(ar n /2G eo )2 

[l + (ar n /2G co )aj 2 + (u*-,^) 


(2.42) 


* * 

The components of the shear modulus G and the compliance J 
for the B. E. L. model are given by: 


4G m (ayn/2G oo ) 3/2 [l + (urn/ZG^)' 2 
{ [l + (ar n /2G oB )2] 2 + (um/ZG^)} 


G* = 


n” 

2G ((BT1/2G ) 

CO ' CO 

r 1 +(a)r)/2G oo ) 2 


'{i 

1 + (o)Tl/2G^) 2 ~\ 

2 + (um/2G : 

CO 

’} 2 


(2.43) 


(2.44) 


J ' = 5 T + 


» ((ur)G / 2)2 


(2.45) 


J" = — + 
(UT] 


(u^G / 2 )‘ 


(2.46) 
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Finally, the dynamic viscosity is given by 


T)* = 


G" 

U) 


n |i 


[1 + ((WI/2GJ2J 


1 .+ (qjr|/2G f _) 2 + (dn/2G )j^ 


/ 


(2.47) 

\/ 


The variations of the storage modulus, the loss modulus and the 
dynamic viscosity with frequency, calculated according to equations 
(2.43), (2.44) and (2.47), are graphically displayed in Figure 2.4. 

As compared with the results for a Maxwell element, displayed in 
Figure 2.2, the B. E. L. liquid model has a longer relaxation time. 

This is consistent with the results of previous correlations based 
upon distributions of Maxwell elements. 

The curves plotted through the data points of Figure 2.3 are 
calculated according to the B. E. L. liquid model from equations (2.41) 
and (2.42). There is excellent agreement with the experimental results. 


2.5 Relationship of Continuous and Oscillatory Shear 

Dyson [ 15 ]] has had considerable success in correlating 

y 

the results of - elastohydrodynamic lubrication experiments with the 
properties of fluids experimentally determined in oscillatory shear. 
Dyson bases his comparison on a simplification of Oldroyd’s [" 39 J 
theory of the steady motion of an idealized liquid. 

The analysis postulates that a simple continuous shear deformation 
includes a rotation of the liquid elements. It is therefore necessary 
to refer all equations that describe its -viscoelastic behavior in 
continuous shear to reference axes which rotate with the element of , 
fluid. The rotating axes yield additional time derivative terms in 
the equation of motion of the fluid and thus additional deformations. 
These equations are solved subject to the velocity boundary conditions. 
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B.E.L. MODEL: Com 
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Figure 2.4. The frequency variation of the modulus components and the dynamic viscosity for the 
B. E. L. liquid model. 






to determine the stresses in the fluid. Finally, the normal stresses 
are described with reference to the fixed axes. 

Dyson's [ 27] simplification of the_ .Pldroy.d-parameters permits 
the normal stresses for a~ f luidTwith^ relaxation time X = r\/G , in 

03 

simple laminar shear, to be expressed in terms of one parameter K: 



2 2 

2_ . KVx 

2 2 2 : 
K 1 + K D X 


P 

.*1 

G 


CO 


1 

K 


KDX 


2 2 2 
1 + KDX 


P = P = 0 


yy zz 

where P = normal stresses 

G = limiting shear modulus 
00 

K = parameter of the analysis 

X = Maxwell relaxation time 

D = shear rate 

x = direction of flow 

y = direction of velocity gradient 

z = direction normal to both x and 

Equations (2.32) and (2.33), repeated below 
for a Maxwell fluid subject to oscillatory shear. 


(2.48) 


(2.49) 

(2.50) 


y 

, have been derived 


£1 = w ^ 2 

G » 1 + m 2 X 2 



00 


ojX 

1 + aA 2 


(2.51) 

(2.52) 


Dyson observed, as a result 
and (2.49) with (2.51) and (2.52) 


of the comparison of equations (2.48) 
, that the shear stress P^ is equal 
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to 1/K of , the va lue of G" at an angular frequency cu = ,KD. Furthermore , 

2 

one half of the normal stress difference, i(P -P ), should be 1/K 

• • ^ xx yy 

of the value of G* at an angular frequency m = KD. A comparison of 
the dynamic viscosity for continuous shear 


* - -t? - 


2 2 2 
1 + K D X 


(2.53) 


with equation (2.34), the dynamic viscosity' in oscillatory; shear, 

G" n 


n* = 


(JD 


2 2 
1 + u> X 


(2.54) 


shows the variation with shear rate D is the same as with angular 
frequency u), with ou replaced by KD. 

The hypothesis above is checked against the results of Russel 
[ 28 ] in Figure 2.5. The variation of apparent viscosity is shown 
for the same three fluids in both oscillatory and continuous shear. 
Note that the two curves begin to diverge at an abscissa value between 
1 and 10. This corresponds to the conditions where G" reaches its 
maximum. 

Whatever model is employed to represent the viscoelastic proper- 
ties of the liquid, its application to continuous shear must be made 
in the rotating coordinate system. Therefore, the generalization of 
this analysis is stated as 


t(D) = (2.55) 

Dyson's application of equation (2.55) to the B. E. L. liquid model 
is compared with the experimental results of Smith [. 29] at low shear 
rates in Figure 2.6. Dyson [ 15..] reports that a constant value of 
K = 7.5 shows good correlation over all Smith’s experimental conditions. 
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Figure 2.5. Comparison of variation of apparent viscosity in oscillatory 
and in continuous shear. 
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Figure 2.6 

Comparison of results of Barlow & Lamb in oscillatory shear with those of Smith in con- 
tinuous shear — mineral oils, steel surfaces. 
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will appear 


As a result of equation (2.55), (G /K) , and not G , 

00 00 

in the equations of motion. A new limiting shear modulus for con- 
tinuous shear is now defined to include the 01droyd-Dyson~ parameter K: ’ 

G = (G /K) (2.56) 

03 CO 

2.6 Limiting Shear Stress 

The Maxwell or B. E. L. liquid model, when applied to continuous 
shear, predicts a shear stress that rises to a maximum and then falls 
with increasing shear rate independent of thermal effects. This be- 
havior is intuitively doubtful and Dyson [ 27 ] reviews the mathematical 
objections. It is suggested that this behavior would give rise to an 
unstable flow pattern. The correlation shown in Figure 2.5 suggests 
a transition to another mechanism of flow as the shear rate approaches 
the value which corresponds to a maximum shear stress. At this shear 
rate, the correlation between experimental and predicted values weakens. 

As an alternative to the falling portion of the shear stress- 
deformation relation, the possibility of a limiting shear stress is 
suggested. The limiting shear stress is the maximum stress a fluid 
can transmit; an increase in the rate of shear can no longer cause an 
increase in the shearing stress. Smith [.29 ] first suggested this 

behavior of a fluid analagous to plastic deformation of a solid. 

Plint's [ 13 ] results further suggest the existence of a limiting 

shear stress in an EHD fluid film. He interpreted the limiting shear 
stress to be the result of a discontinuous shear failure. Dyson [ 15 ] 
suggested the limiting shear stress be a function of the limiting shear 

modulus G . Figure 2.7 shows this results in a good correlation with 
00 

the experimental data of Johnson and Cameron [ 11 3 . 
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Figure 2.7. Dyson's correlation with the experimental data of Johnson 
and Cameron [ 11 J at high shear rates. Curve from 

Dyson [ 15 ] . 
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mean shear stress, MN m" 




Consequently, it is surmised that there are two mechanisms of 
flow for a liquid under the conditions of continuous shear. The 
material properties of the liquid, as well as the transition between 

these t wo mech a nism s/of_j low, are continuous. The liquid jnode 1 for_ 

continuous shear is, therefore, a composite non-linear shear stress- 
strain relation. . It is comprised of a viscoelastic relation for shear 
rates up to the value predicting the maximum shear stress, and a 
limiting shear stress equal to this maximum at higher shear rates. 


2.7 Hyperbolic Liquid Model 

Barlow and Lamb [ 19 ] investigated the viscoelastic relaxation 
in three mineral oils of different viscosity index and composition. 

The experimental results, shown in Figure 2;8, show slight deviations 
from the B. E. L. liquid model. The experimental results are dis- 
placed to higher values of frequency and lower values of the loss 
modulus . 

To add flexibility in the analysis, a "hyperbolic" shear stress- 
strain relation is used which allows easy changes in the limiting 
shear stress or' the rate’of rise to this limit. The relation has the 
additional feature of providing a smooth transition to the flow domi- 
nated by a limiting sheaf stress. The model is mathematically repre- 
sented by 


2 2 



where T = shear stress 

G = limiting shear modulus 


q _ H_ Bii dimensionless shear rate 
G & 

CO 


(2.57) 


(2.58) 
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c = limiting shear stress/limiting shear modulus ratio 
a = rise parameter, rate of rise to limiting shear stress 
decreases as a increases 

Four-mode ls__of. interest, the hyperbolic model for c = .25 and 
c * .20, and the Maxwell and B. E. L. - limiting shear models are 
illustrated in Figure 2.9. 

The hyperbolic liquid model has the following limiting values: 


Lim = 0 

Q - 0 G 


Lim 

Q- 0 







(2.59) 

(2.60) 


Lim = c (2.61) 

fl ® G 


Rbr- ] < 2 -“> 

Q _ oo ■- G J 


For the case of a = 0, equation (2.57) reduces to the true 


hyperbola : 



G 


00 


(2.63) 
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LIQUID MODELS 
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Figure 2.9. Comparison of liquid models. 



CHAPTER III 
TRANSIENT VISCOSITY 



The maximum pressure in the lubricant film between highly loaded 
contacts may be as high as 250,000 psi. The ; lubricant film is there- 
fore subjected to a. large pressure transient as it passes through the. 
contact, and the equilibrium viscosity at the maximum pressure is . 
several orders of magnitude greater than the atmospheric pressure value. 

Measurement of the fractional force between the two contact 
surfaces at low values of slip enables an "effective viscosity" of 
the lubricant to be calculated. At high rolling speeds, this effective 
viscosity is found to be lower than the value calculated from the equi- 
librium value of the viscosity as a function of pressure. The effective 
viscosity also decreases with increasing rolling speed, in a manner 
which is not adequately explained by either viscous heating of ’ the" lubri- 
cant film or by variation of the viscosity as a function of shear rate. 

Fein f 30 ] has suggested the failure of the lubricant viscosity 
to respond to the rapid pressure changes encountered in the contact 
area could be an explanation for the low values of effective viscosity 
which are observed. His analysis shows that under certain conditions 
the time of transit of the lubricant through the contact zone could 
be small compared with the time required for the ' lubricant to reach 
a state of equilibrium following an applied pressure^ step. Consequently, 
the compression of the lubricant never reaches the equilibrium state; 
corresponding to the peak pressure, ahd the viscosity has a lowervalue 
than that measured under equilibrium conditions. An increase in the 
rolling speed reduces the residence time of the lubricant in the contact 
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zone. This results in an even lower value for the : viscosity attained 
by the lubricant, and a consequent decrease in the effective viscosity. 

Chapter III is an analysis of the effect of compressions 1 visco- 
elasticity on the pressure-induced viscosity changes that occur in 
concentrated contact lubrication. The variation of viscosity with 
time, following an applied step in pressure, is described by a non- 
linear model proposed by Kovacs [ 16 ] for the volume creep of polymer 

melts . 

3.1 Compressiona 1 Viscoelasticity 

The response of a liquid to a rapid change in pressure consists 
of an instantaneous volume change, followed by a time-dependent volume 
change. The instantaneous change is attributed to the elastic com- 
pression of the liquid "lattice”, while the time-dependent response 
is attributed to molecular rearrangements. The instantaneous response 
of the liquid, when the experimental time scale is small compared with 
the time, required for molecular rearrangements, is characterized by a 
bulk modulus K . When the experimental time scale is large compared 

CO 

with the molecular rearrangement time, the bulk modulus has a lower 
value, the equilibrium value K^. This behavior may be represented by 
the simple models shown in Figure 3.1. 

Model A is widely used when volume relaxation is investigated 
as a function of frequency. The overall modulus then rises from a low 
frequency. va lue to a high frequency limiting value + K^. 

is the high frequency value of the real part of the complex re- 
laxations 1 modulus, K^(joj) = K^((u) + i K^(m). For model A, the total 
bulk modulus is given by the expression 
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(3.1) 


K = K + K (iuj) = K + K. . — . — 
o r 0 2 1+ 


and the relaxation time 




(3.2) 


where 


r) = volume viscosity 
a) = angular frequency 


Model" B' is more suited to a description of the change in volume 
following a sudden increase in pressure, volume creep, since the 
instantaneous and time-dependent parts of the response are easily 
separated. The response is more simply expressed in terms of the over- 
all compressibility, the reciprocal of the bulk modulus, given by 
equation (3.3) as a function of frequency. 


1 1 


K " T K f (l + iu)X f ) 


(3.3) 


is a modulus associated with molecular rearrangements corresponding 
to changes, in the free volume and is the retardation time given by 


X 

*f K f 


(3.4) 


The viscosity is associated with the changes in the free volume. 
The low frequency or equilibrium modulus is obtained from equation 

i 

(3.3) for cu = 0. 


— + — 
K K £ 

00 f 


(3.5) 


K K £ 

K = ” £ 

o K + K, 


(3.6) 


Models A and B describe the same behavior and a comparison of 


44 



equations (3.1) and (3.3) yields the following additional relations 
between the parameters of the two models: 

K = K + K, 


"t ■ ’’vfe!) 


o 2 
K 2 
2 J 


It follows from equations (3.5) and (3.7) that 


K 


K 


_£ 

K ■ K 


(3.7) 

(3.8) 


(3.9) 


The behavior of liquids is generally found to be more complex 
than that described by these simple models and a combination of several 
models, each with different time constant and moduli, is necessary. 
Alternatively, a continuous distribution of relaxation, or retardation, 
times may be used to characterize the liquid behavior. The introduction 
of a distributed spectrum causes considerable complication in the 
analysis and is not warranted in the present study. Model B of Figure 
3.1 will be used to characterize the behavior of the lubricant. 


3.2 Viscosity Response to a Pressure Step 

The overall change in volume from an initial volume v^ to a 
final volume caused by a pressure change P, is given by the def- 
inition of the secant bulk modulus K : 


V 1 P 


V — V = 

1 2 K 


(3.10) 


The volume change corresponding to the purely elastic deformation 
(v^-v^) is given by 


V 


v -v = 

U i K 


(3.11) 
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Equation (3.10), with the aid of equations (3.5), 
may be written as 


(v 


rV 


+ 


<v.-v 2 ) - 




(3.10) and (3.11), 


(3.12) 


Therefore , 


v . 
i 


v 




(3.13) 


Equation (3.13) may be taken as a definition of K^. 

The time dependence of this volume change is given by the 
parallel spring and dashpot combination of model B. The response is 
governed by 


P 


dv 

dt 


+ K, 



Combining equations (3.13) and (3.14) yields 


(3.14) 


^f dv 


K f dt V 2' V 


(3.15) 


where v varies between and For small changes in pressure, when 

and can be regarded as constants, equation (3.15) has the solution 
v-v 2 = ( v i“ v 2) expC-t /\ f ) (3.16) 

and the total response to the pressure step is 


7~ = + - ex P (-t/X f )] } 

1 00 £ 


V v 

v 


(3.17) 


where the retardation time 

However, for large pressure changes the parameters ru , K and 
can no longer be regarded as constants. In particular, the vis- 
cosity r|j may be expected to change by many orders of magnitude under 
the pressures occuring in the contact zone. There is considerable 
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evidence from ultrasonic studies of liquids that the volume viscosity 
T] , and hence , is closely related to the shear viscosity r) and has 
the same temperature dependence. Litovitz and Davis [ 3] ] and 
Taskoprulli, Barlow and Lamb f 45 j offer such-evidence for Liquids 
including lubricating oils. It is assumed here that has the same 
dependence on the free volume as the shear viscosity. Then is 
related to the free volume v^ by the Doolittle [ 24 ] equation: 

In n f = A + B/f (3.18) 


v-v 

where f , fractional free volume (3.19) 

o 

v = specific occupied volume 
A,B = constants 

The value of A is characteristic of the liquid; the value of B is 

usually close to unity. The occupied volume, a function of pressure, 

is assumed to be independent of time and thus is associated with the 

instantaneous bulk modulus K . The variation of with pressure is 

00 t 

described by the parameter s, defined by 


s = 



(3.20) 


where f is the final volume of the fractional free volume, and r) f 
l t 2 

is the final value of r^, at pressure P. If no change in the fractional 
free volume occurs during the instantaneous compression, the initial 
va lue of s is 


s, = 



(3.21) 
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where n f and f are the values of n f and f at the initial equilibrium 
f 1 1 t 

state when P = 0. The shear viscosity is also described by equations 
(3.18) and (3.20), so that. 


■i - 


(3.22) 


where and are the initial and final equilibrium values, respectively, 
of the shear viscosity n • 

Following Kovacs [ 16.] , if it is assumed that the occupied 
volume v^ remains constant after the initial compression, the parameter 
s is evaluated from its definition: . : 


„ r (V ' V2)V ° i 

~ L(v-V )(v -V )J 


(3.23) 


o £ o 

Therefore, the differential dv is 

(v-v ) 2 

dv = ds 


(3.24) 


Equation (3.15) is written in terms of the parameter s with substi- 
tutions from equations (3.20) and (3.24): 


_ exp(-s) (v -V . dt 
s (v 9 -v ) “ x 9 


(3.25) 


Noting that 


1 - 


Hi = iV2 ' V o ) 

B . (v-v ) 
o 


equation (3.25) becomes 

exp(-s) 1_ 




. dt 

ds = - — 


^2 


(3.26) 


(3.27) 


where i- s a retardation time characteristic of the final equilibrium 
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state , given by 



ex P( s 1 > 


(3.28) 


The term (sf 2 /B) in equation (3.27)_is typically much less than-nriity, 
so that the expression (l-sf^/B) ^ may be expanded to give 


. ex gCr s ) ds + exp(-s) ~ ds = - 
s r B -> . 


dt 


(3.29) 


as the differential equation describing the time-dependent compression 
of the liquid. This may then be integrated from the initial value s^ 
at t = 0 to an intermediate value s at time t: 


Ei(-s 1 ) - Ei(-s) + ^ |^exp(-s)-exp(-s 1 )J = 


0.30) 


where Ei is the exponential integral. The viscosity at time t is 
given by 


•q = n 2 exp(-s) 


0.31) 


Figure 3.2 shows the variation of viscosity with time predicted 
by equation (3.30) for the following parameters: P = 200,000 psi; 

rij = 10 ^lbf-sec/in^; n 2 = 10^1bf-sec/in^ ; f^ = 0.05; B = 1. These 
values are typical of those experienced by a lubricant in the contact 
zone of a heavily loaded rolling contact. For values of t of the order 
of X 2 or less, the viscosity is seen to be significantly less than the 
equilibrium value. 


3.3 Viscosity Response of a Lubricant to a Pressure Step 

To determine whether the behavior described by equation (3.30) 
has a significant effect on the properties of the lubricant, values 
of the residence time of the lubricant in the contact and the retardation 
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Viscosity (lb f sec in ' 2 xIO 2 ) 




time must be determined. The lubricant in a rolling contact is 
subjected to high pressure for a time equal to 2b /U where b is the 

half-width of the contact zone and U is the rolling speed. For typical 

-2 — - 

values of b = 10 in and U = 100 in/sec, the residence time is of 

-4 

the order of 10 sec. 

The time constant ^ is characteristic of the final equilibrium 

state of the liquid. Values of and for lubricants are not 

available, but reasonable estimates may be made from ultrasonic data 

on other liquids. The viscosity ri£ is related to the volume viscosity 

2 

Tl v by equation (3.8), n f = ^(K^/K^) . Litovitz and Davis [ 31 ] 

report that the ratio K /K 0 is of the order of 3 for many liquids. 

co Z 

2 

Therefore, a value for (K /K ) of 10 may be used. Ultrasonic studies 

co Z 

also indicate that n v is closely related to the shear viscosity; a 

ratio of u v /ri = 5 has recently been reported by Barlow, Lamb and 

Taskoprulu [ 32 ] . The value of is then given by 50^, where 

is the atmospheric pressure shear viscosity. The time constant i s 

given by expCspn^/K^ = SOr^/K^. The bulk modulus is related to 

the relaxational modulus by equations (3.7) and (3.9); for a ratio 

K f /K Q = K^/I^ = 3, then = 61^ . But is experimentally found to 

be approximately equal to 4/3G , where G is the high frequency limiting 

00 00 

4 

shear modulus of the liquid. A value for G of 4.35 x 10 psi 

00 

9 -2 

(3 x 10 dyn. cm ) has been reported by Hutton f 33 ] for a H.V.I. 
lubricating oil at 30 °C giving a value of of 3.5 x 10"* psi. This 
modulus will change significantly with pressure. Dyson f 15 J reports 
measurements of G as a function of pressure give a typical value for 

CO 

$G /9P of 3, and K for lubricants varies in a similar manner in the 
oo o 

Pressure-Viscosity Report [ 34 ] . If it is assumed that the ratio 
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K^/K q remains independent of pressure, then = (3.5 x 10^ + 9P) psi, 
and the retardation time is given by 

. son, 

\ 2 £ t (3.32) 

(3.5 x 10 + 9P) , . 

For the values given above, \ 2 has a value of the order of 
_2 

2 x 10 sec, which is much greater, than the residence time of the 
lubricant in the contact zone. The "instantaneous viscosity" of the 
lubricant will therefore be much less than the equilibrium value, 
resulting in greatly reduced values of effective viscosity in accordance 
with experimental observations. 


52 


CHAPTER IV 


MATHEMATICAL FORMULATION 


The present analysis of traction in e la stohydrodynamic contacts 
includes the iterative solution of the momentum and energy equations 
with the fluid properties functions of pressure and temperature. The 
shear rate and transient time effects have been isolated as discussed 
in Chapters II and III. 

In this chapter the momentum and. energy equations are developed 
and the pressure profile, the film thickness and the material property 
functions are discussed. The set of equations developed are then 
solved numerically. - - 

4.1 Geometry and Coordinates 

..... 

The geometry of a typical disk machine is ; shown in Figure 4.1. 

Two cylinders of radii and , rolling with velocities and U^, 
respectively, are separated by a lubricant film of thickness 2h. A 
closer view of the contact zone as shown in Figure 4.2 is more useful 
for the purposes of this analysis. The disks have deformed elastically 
to form a contact zone of width 2b and the film thickness is approxi- 
mately constant with the surfaces of the disks remaining nearly parallel. 

The coordinate system is defined to have the origin on the center 
lines of both the fluid film and the flat contact zone. The x-axis 
is the center line of the lubricant film with the positive direction 
in the direction of flow; while the y-axis, the perpendicular bisector 
of the flat contact zone,. is arbitrarily taken positive toward the 
disk rolling with velocity U^. The z-axis, not shown in Figure 4.2, 
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Figure 4.1. Typical disk machine geometry. 
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is perpendicular to both the x and y-axes with the positive direction 
consistent with a right-handed Cartesian coordinate system. 

The control volume of interest is defined as an element of fluid 
of length dx in the direction of flow, bounded by the disk surfaces 
in the y-direction and of unit thickness in the third direction. 


4.2 Pressure Distribution 

The pressure distribution in the contact zone is assumed to 
have the elliptical Hertzian dry contact profile given by 

- %z 7 1 - (!) 2 < 4 -‘> 

where p = maximum Hertzian pressure 

x = distance from the center of the contact 

4R P Hz - 

b (4.2) 

5 

= half Hertzian width 
R = effective radius of the disks 
E = effective modulus of elasticity 
The deviations from this assumed distribution are mainly in the entrance 
zone at low pressure levels. Their effect on the sliding friction is 
very small and is neglected. 


4.3 Film Thickness 

The minimum film thickness in elastohydrodyriamic contacts at 
moderate rolling speeds can be accurately predicted by the Dowson 
and Higginson [ 35] formula: 


h 

o 


1.6 a 0 ' 6 (p D )°' 7 ( E ) 0 - 03 R 0 ' 43 
ent 


w 


0.13 


(4.3) 
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where h = minimum film thickness 

o ‘ ■■ 

a = viscosity-pressure exponent 

n = viscosity of the lubricant at the conditions of 
'ent - - -■ 

entry to the contact 
U = mean rolling speed 

E = effective modulus of elasticity 

R = effective radius of the disk pair 
w = load per unit length of cylinder 
Note that the minimum film thickness is only slightly dependent on 
the load w and virtually independent of the elastic modulus E. 

Dowson and Higginson [ 36] suggest the parallel film thickness 

2h is 2CK greater than the minimum film thickness h^. 

The Dowson and Higginson prediction of film thickness is based 
on an isotherma f ana lysis which is no longer adequate for heavily loaded 
contacts operating at high rolling speeds. Cheng [ 37 ] has calculated 
the lubricant film thickness in the Hertzian flat for high speed and 
heavily loaded rolling and sliding contacts. He used a Grubin-type 
inlet analysis including full thermal-hydrodynamic effects. The results 
obtained for a wide range of load, speeds and lubricant properties 
showed that the loss of film due to thermal effects is strongly in- 
fluenced by the rolling velocity and the inlet viscosity of the lubri- 
cant, while it is somewhat insensitive to the change of load. The 
presence of sliding does not have a significant influence on- the calcu- 
lated film thickness, whereas the rolling speed has a far more pre- 
dominant effect at the inlet. -v • 

The loss. of. film thickness due to therma 1 effects can be most 
conveniently represented by a thermal reduction factor §^,, which is ’ 
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defined as the ratio of the actual film thickness to that predicted 
by isothermal theory. Run #29 of Cheng’s work is most applicable to 
the lubricant properties of this study and has been reproduced in 
Figure 4.3. Cheng's parameter is defined by 


where 


2ri . u 

Q .Jlsat — 

to k I 

ent 

2 

•n = viscosity at entry conditions (lbf-sec/in ) 
'ent 


(4.4) 


U = average rolling speed (in/sec) 

k = thermal conductivity of lubricant (Btu/°F-hr-f t) 

T = temperature of lubricant at entry (°R) 

The film thickness including the thermal effects is calculated 
by multiplying the isothermal film thickness, based on the Dowson- 
Higginson formula, by the parameter determined in Figure 4.3. For 
example, at a rolling speed of 500 in/sec at 175 °F, equation (4.4) 
requires 

n 2(.87 x 10" 5 )(500) 2 _ „ 

m ( .0216) (635) 

and Figure 4.3 determines the thermal reduction factor 


«! = - 8 1 

Values for other conditions are similarly calculated. The results 
are shown in Table 4.1. 


4.4 Momentum Equation 

In applying the principle of conservation of momentum to the 
lubricant in the contact zone, we make the following assumptions: 
1. For the case of a line contact, all the variables are 
independent of z, the direction of the axes of the disks. 
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Table 4.1 

VALUES OF THERMAL REDUCTION FACTOR AT 
TYPICAL EXPERIMENTAL CONDITIONS 



175 °F 

220 °F 


* 32 

r-H 

• 

II 

cr S 

500 in/sec 

= .81 

$ T = .89 


Q = 1.27 
m 

= * 60 

1000 in/sec 

§ T = .58 

$ T = .72 








2. As compared with the lubricant film thickness, the radii of 
curvature of bearing components are generally very large. In the 
specific case of the disk machine, the radii of the disks and 

» 2h. Accordingly ,“a"lleff ects due to curvature of the fluid film 
are neglected. 

3. : As compared with the much larger pressure and viscous forces, 
the inertia and body forces of the lubricant are negligible. This 
imples that the pressure and viscous forces acting on the fluid are 

in equilibrium. 

4. As compared with the other dimensions of a lubricated con- 
centrated contact, the film thickness is very small. Therefore, the'.-' 
derivative of u with respect to y is large in comparison with all 
other velocity gradients. 

5. The pressure gradient across the lubricant film is also 

insignificant due to the relative smallness of the film thickness. 
Accordingly, p = p(x) ± p(x,y). ; 

The assumptions outlined above reduce the surface forces acting 
on a fluid element in the contact zone to those shown in Figure 4.4. 

The momentum equation can be derived directly from the balance of 
these surface forces. Equilibrium in the x-direction requires 


J*E = 
dx 



(4.5) 


The shear stress j must have two components. The rolling of 
xy 

the two cylinders produces the first component, while the second 
component results from the difference in. rolling velocities, or slip. 
Consider a control volume bounded in the y-direction by the surfaces 
of the disks as shown below. 
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V p + 6x ; 6y 



Figure 4.5. Forces acting on the control volume. 




The size of the fluid element considered in Figure 4.4 can be 


increased to that of the control volume used in Figure 4.5 by inte- 
grating equation (4.5) with respect to y over the film thickness. 


h 


‘ -h 


dy = f + h dy 
dx - _ h ay 


or 


2h^ 

dx 


dx 


= (t v „) - (t 


Xy y=h 


Xy y= -h 


“ T 


roll 


(4.6) 


Thus, the rolling component is independent of the slip and is a function 
of the pressure gradient through the contact zone. Only the component 
of stress that arises due to the relative sliding of the two disks is 
of interest in this study. This component is easily separated by 
neglecting the pressure gradient term of equation (4.5). For convenience, 
we redefine ( T ) = t and the final form of the momentum equation 


'xy 


slip 


becomes 


ZL = 0 

dy 


(4.7) 


The shear stress j is supplied by one of the rheological models 
considered for the lubricant. The form of the model is 



Therefore at any position x 


(4.8) 


and 



(4.9) 

(4.10) 


where is a constant of integration. When a specific model is used, 
one can isolate ^ and integrate with respect to y from y = -h to y = h. 
Moreover, it is assumed that the profile for ^ is symmetric with respect 
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Figure 4.6. Energy, balance. for a fluid element. 
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It is assumed that the heat generated in the fluid element will 
be dissipated in two modes. Convection, the first mode of heat transfer, 
is the process by which energy is carried out of the contact zone with 
the lubricant. In the second mode, heat will be conducted across the 
film to .the disks. Since the lubricant film thickness is small in 
comparison with the Hertzian contact width (x-direction) and even smaller 
in comparison with the cylinder width (z-direction) , the temperature 
gradients in the x and z-directions must be small in comparison with 
those across the film. Therefore, only conduction in the y-direction 
is considered. These modes of energy transfer are shbwn for the fluid 
element in Figure 4.6. 

The rate of heat generation per unit volume q is given by the 
product of stress with the rate of strain. 

q = T 2H (4.13) 

T By 

An energy balance on this fluid element requires that the net energy 
into the control volume be zero. 

- p c - p c + k - T = 0 (4.14) 

p s* p sy sy 2 sy 

Heat Transported to Heat Conducted to Heat Generated in 

Control Volume Control Volume Control Volume 

where T = temperature of the lubricant 

p = density of the lubricant 

c = specific heat of the lubricant 

k = thermal conductivity of the lubricant 

Continuity requires that 

^ + ^ = 0 (4.13) 

5 X By 
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Therefore , 


- ,cL& + + k^| - - T aa (4.i6) 

p V a x ay 1 2 ay 

Convection Cross- Conduction Viscous Heat Generation 

Convection 

The ratio of convection to conduction is estimated, by assuming 
a triangular temperature profile, to be 


pcUh 2 

2bk 


(4.17) 


Equation (4.17) demonstrates that convection will have its 

2 

largest effect for a maximum value of (Uh /b). This corresponds to the 
condition of maximum rolling speed and minimum load. 

For the thin lubricant films in EHD contacts, where h« b, the 
convective heat transfer can usually be neglected. Therefore, the 
governing energy equation may be written as a balance of viscous heat 
generation and heat transported by conduction. 




(4.18) 


The consequences of this assumption are discussed in section 6,3. 

The lubricant in contact with the disks assumes the surface 
temperature of the disks. Blok [ 38 ] has analyzed the problem of a 
moving heat source. His results demonstrate that the disk surfaces 
in the concentrated contact will have a mean "flash temperature" 
higher than the bulk temperature of the disk. Equation (4.19) is 
the expression Blok derived for the flash temperature. 

0.48 nw|U r U 2 | 


T -T, = 
s b 


(k p c U b) : 
nrm m 


(4.19) 
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where 


T g = mean surface temperature in the contact zone, 
"flash temperature" 

T, = surface temperature entering the junction, bulk 
b 

temperature of the disk 
k^ = thermal conductivity of the disks 


c = specific heat of the disks 
m 


n = density of the disks 
r m 


4.6 Equilibrium Viscosity Function 

The equilibrium viscosity is the viscosity measured after the 

lubricant has reached a state of static equilibrium under a given 

temperature and pressure. Viscosity deserves special attention in 

the study of friction in concentrated contacts. Unlike other physical 

properties, which change only slightly with temperature and pressure, 

the viscosity of a lubricant can change by several orders of magnitude. 

Viscosity is most simply defined by Newton's Law: 

t = T|V (4.20) 

2 

where X = shear stress (dynes /cm ) 

y = shear rate (sec '*’) 

T) = viscosity (Poise) 

This can be generalized to equation (4.21) for a viscoelastic fluid. 



(4.21) 


The viscosity of a liquid is basically the resistance of molecules 
to move past the force fields of neighboring molecules. It is a compli- 
cated pressure and temperature-dependent function. 

The viscosity of a liquid and the rate of change of the viscosity 
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due to a temperature change decrease with increasing temperature. An 
increase in temperature of the fluid increases the thermal agitation 
of the molecules which, in turn, lessens the forces of attraction be- 
tween molecules. Thus the viscosity decreases. 

The has been considerable effort to find an accurate relationship 
for predicting the variation of viscosity with temperature. Some of 
these relationships have theoretical foundations but the empirical 
formulas provide the most satisfactory predictions. The viscosity- 
temperature relationship found by Herschel [ 46 ] is 

to8l0 fe) ' P ' 1C8l ° T <4 ’ 22) 

where ri is the viscosity (centipoise) at the temperature T (°F) and 
n ref and P are constants. Thus the "Herschel Chart", a plot of equa- 
tion (4.22) on log-log graph paper, is a straight line for a given 
lubricant. The equation is simple but Appeldoorn [ 40] has found it 

surprisingly accurate for oils of very different viscosities. 

The viscosity-temperature data for Mobil XRM 109 F4 and Shell 
Turbo 33 is given in Table 4.2. This data determines the Herschel 
equations : 

log 10 n = 8.974 - 3.2 log 10 T (4.23) 

for Mobil XRM 109 F4 and 

log 10 n = 7.3409 - 2.8 log 1() T (4.24) 

for Shell Turbo 33. Both of the equations above may be plotted on 
log-log graph paper as shown in Figures 4.7 and 4.8 and used as con- 
venient Herschel Charts. 

The effect of pressure on viscosity is influenced by both the 
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TABLE 4.2 


VISCOSITY -TEMPERATURE DATA 








Figure 4.7. Herschel Chart for Mobil XRM 109 F4 calculated from 
equation (4.24). 
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Figure 4.8. Herschel Chart for Shell Turbo 33 calculated from 
equation (4.25). 
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pressure level and the bulk viscosity of the fluid. The same increase 
in pressure will have a greater effect on the viscosity at a high 
pressure level than at a lower level. This results from the fact that 
more of the free space between molecules is already taken up at the 
higher pressure level. The same effect is responsible for a fluid of 
high viscosity undergoing a greater viscosity change than a fluid of 
lower, viscosity for the same increase in pressure. 

The Pressure-Viscosity Report [ 34 1 , which includes data 
on several paraffinic and napthenic mineral oils, pure hydrocarbons 
’ and synthetics, is an excellent source of pressure-viscosity data. 

Chu and Cameron [ 41 ] have analyzed the results of this report 
in an attempt to find a sufficiently accurate pressure law and corre- 
lation. The usual simple exponential law was found inadequate for 
paraffinic oils. Paraffinics were found to obey the law 

(log 1Q ri ) 3 ^ 2 = m(p + a) (4.25) 

and there was a simple correlation between m and n. the base 

'base 

viscosity. Including this correlation, equation (4.25) becomes 

2/3 

log 10 n = 0.18(log 10 r ^ se ) 2/3 ( P + 13.2 Jlo &lQ n base ) (4.26) 

where r) = viscosity in centipoise at pressure p 

tv = base viscosity at p = 0 
'base r 

- p = pressure in kpsi • - 

Note that this convenient form of the Chu and Cameron viscosity-pressure 
law automatically correlates to each lubricant through . 

Cheng [ 37] has analyzed, the data of the same. PressurerViscosity 

Report. He used the following alternative viscosity-pressure relationship: 
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■ exp^ap + (P + YP>(f - ] 


(4.27) 


where r) = viscosity at pressure p ‘ ~ _ 

r| = reference to viscosity at p = 0 and T = 
p = pressure (psi) 

T = absolute temperature (°R) 
a = viscosity-pressure coefficient 
P =,.5.1 x 10 7 a 

V = 930 a" - 

Figure 4.9 exemplifies the pressure dependence of the equilibrium 
viscosity function according to. Cheng, and Chu and Cameron. It has 
been calculated for the Mobil XRM 109 F4' lubricant at 175 °F. 


4.7 Limiting Shear Modulus 

The pressure and temperature function for the high frequency 
limiting shear modulus has been developed by Dyson [ 15 ] in a corre- 
lation with Smith's [ 29 ] experimental data. The development is 
outlined below. ’ 

l ~‘\ 

Hutton [ 33 ] experimentally determined that a high viscosity ; 
index mineral oil at atmospheric pressure varies with temperature ac- 
cording to , . .. 


— = 2.52 + 0.024 T 
00 


(4.28) 


and referenced to conditions at 20 °C, 

G » (T> 3 " • • ‘ ‘ 

G (20°C) 2,52 + 0,024 T 

CO 

where G is in GNm 2 (10^° dynes /cm^) and T is in °C. 

CO 


(4.29) 


73 



Log lo ij(cp) 



equations (4.26) and (4.27). 
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Variation with pressure of the high frequency limiting shear 

modulus is more difficult to estimate, since there is insufficient 

lubricant~da ta — at~h'i'gh pressures . Guided by equation (4.29), Dyson 

looked for a correlation of the shear modulus with the quantity 

3 p 

2.52 + 0.024 T 

Figure 4.10 is the correlation found with Smith’s experimental results. 
Although this correlation predicts an impossible negative value for G 

03 

at p = 0, the predicted values at higher pressures are the best avail- 
able. The limiting shear modulus function, as determined from Figure 
4.10, is 

G co^P ,T > = °* 4 [ 2.52 + P 0.024 t] " 10 (4.30) 

Converting equation (4.30) into English units, one obtains 

G >’ T) = 2.52 + *.bl33(T-492) “ 1,45 x 10 (4.31) 

where G^ is now in psi, p is in psi and T is in °R. Equation (4.31) 
has been used in determining the limiting shear stress in the liquid 
mode Is . 

4.8 Numerical Solution 

The numerical solution of the equations governing the friction 
in elastohydrodynamic lubrication is a Fortran IV coding for use on 
a CDC 6400 digital computer. The program is outlined in Figure 4.11 
and a complete listing is given in Appendix B. 

Program CONTROL is the backbone of the traction calculation 
calling on several subroutines as they are needed. The load is a 
Hertzian elliptical pressure profile, as developed in section 4.2, and 
the lubricant film of uniform thickness is calculated according to 
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: 2.-52 + 0-024T 


Figure 4.10. Correlation of limiting shear modulus with- the experi- 
mental values of Smith [ 29 ]. Curve from Dyson [ 15 ] . 
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VISC, RTMI, ZERO, EXPI 


i i 





Program CONTROL and its subroutines. 























section 4.3. The material properties of the lubricant are allowed to 
vary as functions of local pressure and temperature as discussed in 
sections 4.4 and 4.5. The variation of the lubricant properties also 
includes the shear rate effects proposed in Chapter II and the transient 
time dependence analyzed in Chapter III. An iterative solution of the 
momentum and energy balances (sections 4.6 and 4.7) is used to determine 
the shear stress, and the velocity and temperature profiles in the 
contact zone. The tractive force on a disk surface, resulting from 
a given sliding velocity, is determined by integrating the shear stress 
at the disk surface over the contact area. The traction coefficient 
is then defined as the tractive force divided by the applied normal 
load. 

Function VISC supplies the transient value of the viscosity 
according to the model analyzed in Chapter III. For the purposes of 
computation, the contact area is divided into six equal zones, the 
pressure being taken as constant within each zone. For the first two 
zones, the pressure step is assumed to be applied at the beginning of 
each zone. The viscosity attained at a time corresponding to passage 
through half the zone is calculated, and this value is used as an 
average viscosity for the zone. For a rolling speed U and contact 
width 2b, this time is b/6U sec. For the third zone, allowance is 
made for the viscosity increase in the preceding zones by calculating 
the viscosity at a time b/2U sec. In each case, the initial viscosity, 
at the instant of applying the pressure step, is taken as the viscosity 
at atmospheric pressure and the disk temperature. For simplicity, a 
viscosity distribution which is symmetrical about the center of the 
contact is assumed, although the actual distribution is asymmetrical. 
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with the maximum viscosity occurring on the exit side of the center. 
This approximate method provides a rapid and simple method of computing 
the effect of the rolling speed on the viscosity of the lubricant. 

If the transient effect is to be neglected for any reason, the 
following trivial subprogram may be substituted for Subroutines VISC, 


R1MI, ZERO and EXPI. 


FUNCTION VISC (P ,ETA2,C0DE) 

VISC = ETA2 

RETURN 

END 


Subroutine RTMI supplies the solution s^ of equation (3.30) 


EiC-s^ 


t 2 r 

- Ei(-s) + — exp(-s) 



to Function VISC. Muller’s iteration scheme of successive bisections 
and inverse parabolic interpolations is used. RTMI is available in 
the Vogelback Computing Center at Northwestern University. Its listing 
is included in Appendix B for completeness. It is a requirement of 
RTMI that equation (3.30) be represented as a separate function sub- 
program. Function ZERO meets this need. 

The evaluation of the exponential integral in equation (3.30) 
is performed in Function EXPI. This routine computes the exponential 
integral for negative arguments in the range -20 to zero. For negative 
values of argument x the exponential integral is defined by 


Ei(x) = J 

J -x 



(4.32) 


In equation (4.32), a polynomial approximation is obtained, for values 
of the argument between zero and -5, by means of the Taylor series 
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expansion by Luke and Wimp [ 42 ] ? 


where 


14 

EXPI(x) = ln|x| - £ b v (-x) V 

v=0 



= - .57721566 
= 1.0 
= - .25 

= .055555520 

= - .010216662 
= .0016666906 


= - 

.23148392 

X 

10- 3 

= . 

.28337590 

X 

10" 4 

= - 

.30996040 

X 

10 -5 

= 

.30726221 

X 

10‘ 6 

= - 

.27635830 

X 

10" 7 


.21915699 

X 

IQ" 8 

= - 

.16826592 

X 

10" 9 

= 

.15798675 

X 

io - 10 


.10317602 

X 

10- 11 


(4.33) 


Equation (4.34) is the exponential approximation used for . arguments 
in the range -5 to -20. 

■ ' ‘ x " 

EXPI(x) = -2.658760 -3 (4.34) 

Function PSI specifies the shear stress^strain relationship to ; 
be used in the momentum equation. The liquid model may be changed 
simply by replacing the deck of this function. Routines for the fol- 
lowing three liquid models are inc luded in the listing of > the . program:.; 
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' ; 1. Maxwell- - Limiting Shear Stress Model 

2. ^ Bv E. L. - Limiting Shear Stress Model 

3. Hyperbolic Shear Stress-Strain Model. 

Hie routine defines the function "Y as . 


. 3 


2 l 


Y = log 


10 


h M' d y 

o a? 

<VV 


(4135) 


where the velocity gradient gu/gy is a .function of the shear stress, 
and therefore, dependent upon the liquid model. , 


_ -i -./Wf-v-- 


gu °° 

Sy ri " * 


* > 


2 3- 


G 

00 


for the Maxwell model; 


(4.36) 


~ — — - f 5 5 . 2 (, \ +4-1 • 

ay n L V a ) a J 


(4.37) 


for the B. E. L. model; and 


c.-- 


_ G G 

gu _. «> , °° 

n 


(4.38) 


for the hyperbolic model.' The boundary conditions specified by equa- 
tions (4.11) and (4.12) require 


f ^dy = 

• J'n' 3y 


(u 2 - u i) 

: 2 ,■ 


Equations (4.35) and (4.39) combine to require 

Y =0 


(4.39) 


(4/40) 
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Subroutine SECANT determines the shear stress solution of the 


momentum equation by solving equation (4.40). This routine is a modi- 
fication of Newton's method. For any function f(x), two initial guesses 
of the root x^ and x^ are required. A straight line is "drawn" through 
f(x^) and f(x£) and extended to cross the x-axis. This new point x^ 
determines f(x^) which is then connected with fCx^) to determine x^, 
as depicted in Figure 4.12. This process continues until the root of 
the function is determined. A number of checks are also included in 
Subroutine SECANT to both guarantee and expedite convergence. 

Subroutine INTEG integrates any non-equidistantly tabulated 
function f(x^) between the limits a and b, where a or b must equal 
f (x^) . The integrated function may be defined as 

INTEG ff(x.)l = [ f(x.)dx. 

L a. J x 1 

A method of overlapping parabolas is employed with suitable modifications 
to yield the fastest possible integration with second order accuracy. 

The development of the quadrature for this subroutine is shown in 
Appendix A. Subroutine INTEG is called upon to integrate the velocity 
gradient in the solution of the momentum equation and again, to inte- 
grate both the Laplacian and temperature gradient in the solution of 
the energy equation. 

Subroutine PRINTS provides the numerical output of the program, 
and Subroutine SETPLOT issues rapid line printer plotting of the 
temperature profiles and traction coefficient versus slip curves. SETPLOT 
is a library routine of Vogelback Computing Center and will probably 
require considerable changes in the coding for use at another facility. 
Since it is not required for the traction analysis, the listing for 
SETPLOT is not included with that of the program. 
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CHAPTER V 


EXPERIMENTAL INVESTIGATION 


The design and manufacture of a disk machine was completed as 
part of this study. The purpose of the experimental investigation 
was to gather extensive data for two new synthetic lubricants, 

Mobil XRM 109 F4 and Mobil XRM 177 F4. The conditions under study 
were those of high loads and high rolling speeds where there was. a 
paucity of experimental data. Special emphasis was given to the effect 
of additives upon the frictional torque. This chapter describes: the 
disk machine, the lubricant properties and the test procedure. 

5.1 The Disk Machine 

The design of the disk machine for this experimental investigation 
was guided -by the following requirements. The disk machine must be 
capable of accurately measuring the tractive force transmitted across 
the line contact of the two disks for a wide range of loads, rolling 

: t 

speeds and slips. A sufficient normal load is required between the 
disks to insure operation in the e la stohy dr odynamic . regime 1 The drive 
to the disks must allow easy adjustments of the mean rolling speed 
and the amount of sliding at the contact. The lubricant must be de- 
livered to the contact at a controlled rate and temperature, instru- 
mentation is required to measure the normal and tangential forces on 
the disks. The angular velocities of the disks, as well as the slip 
or difference in the disk velocities, must also be accurately measured. 
Finally, the surface temperature of the disks as they enter the con- 
tact zone is required for an accurate knowledge of the friction. A 
detailed description of the machine designed to meet these requirements 
follows . 
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The disk machine, pictured in Figures 5.1 through 5.3, was 
designed with two 6-inch diameter disks. These large disks were se- 
lected to allow a direct drive system at high speeds, thus minimizing 
any possible vibrations. The lower disk is supported on two high- 
speed roller bearings which are mounted in the main frame and the 
upper disk is contained in a loading arm which is hinged on the frame 
with a spherical roller bearing. 

The load is applied by an air cylinder at the far end of the 
loading arm. With a 30 psi air supply, a maximum Hertzian stress of 
300,000 psi can be obtained for a £-inch contacting width. The ap- 
plied normal load is monitored by a four-strain gauge bridge mounted 
on the air cylinder shaft. This is necessary for accurate measure- 
ment of the normal load, as the friction in the air cylinder is incon- 
sistent. The loading arm was designed to permit the necessary align- 
ment to insure a uniform load across the line contact in the axial 
direction. 

Each of the disks is attached through flexible couplings to 
separate 40 hp D.C. field controlled electrical machines. The shunt- 
field current method of speed control is simple and efficient and the 
speed regulation, for a given speed adjustment, is excellent. The 
complete electrical circuit, schematically shown in Figure 5.4, is 
the Hopkinson mechanical-loss-supply feedback circuit described, for 
example, by Kloeffler, Kerchner and Brenneman [ 43 ^ . This circuit 
was inspired by a related feedback system successfully used in the 
experiments of Jefferis and Johnson [ 12 ] . One of the D. C. 

machines behaves as a motor driving one of the disks. The second disk 
is driven by the friction force transmitted at the contact and drives 
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Figure 5.3. The disk machine, 
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Figure 5.2. The disk machine 
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ARMATURES CONNECTED BACK-TO-BACK 
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Figure 5.4. Electrical connections for the disk machine 







the second D. C. machine as a generator. The energy is electrically 
recycled by supplying the generated voltage to the armature of the 
motoring machine. The losses in this cycle are mechanically supplied 
by a 20 hp motor connected to the double-ended armature shaft of the 
40 hp D. C. motor. 

The speed of the upper shaft is controlled by the small 20 hp 
booster motor. A wide range of rolling speeds can be obtained by 
varying the supply voltage and the field resistance of the booster 
motor. The speed differential between the two shafts is controlled 
by adjusting the field resistors of the 40 hp D.C. motors which are 
connected together across the armature terminals. This arrangement 
allows a wide and continuous range of sliding speeds to be easily 
obtained. This differs greatly from the typical two-disk machine 
arrangement in which the slide-to-roll ratio is fixed by a gear ratio, 
or the ratio of the diameters to the disks. 

Jefferis and Johnson reported torsional vibration difficulties 
with the original design of their disk machine. Every effort was made 
to eliminate vibration problems in this design. The lower shaft of 
the disk machine can be approximated as indicated in Figure 5.5, for 
the purpose of determining the natural torsional frequency of the shaft. 
The lowest natural frequency of this system is calculated to be 55 
cycles per second which is slightly above the condition existing at 
the maximum experimental rolling speed of 1000 in/sec. The flexible 
couplings used in the system were chosen, in part, for their high 
damping characteristics to further insure smooth operation. Vibration 
problems were not encountered in the course of the experiments. 

The frictional torque transmitted through the line contact is 
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SIMPLIFIED MODEL 







measured by a torquemeter mounted in the lower shaft of the disk 
machine. The measured torque and the known radius of the disk are 
used to calculate the tangential tractive force of the contact. The 
torquemeter consists of a four-strain gauge bridge mounted on a cali- 
brated torsion shaft. The electrical output signal of the torquemeter, 
passed from the rotating shaft through a set of slip rings, is a measure 
of the instantaneous torque in the lower shaft. This signal is dis- 
played on an oscilloscope or measured by a digital voltmeter if an 
integrated average value is desired. The torque measured in this 
manner includes the frictional torque of the two lower support bearings 
which is accounted for as follows. The Hopkinson electrical circuit 
allows rotation of the disks in both directions and, therefore, pure 
rolling is possible. The bearing friction and any minute rolling 
friction are calculated from the measurements at pure rolling. The 
combined value is small and averages to a frictional torque corre- 
sponding to a friction coefficient of 0.002. 

The angular velocity of each disk is measured by a timing wheel, 
seen in Figure 5.3. Each timing wheel has 100 equally spaced holes 
along the circumference. Light supplied from a high intensity source 
to the timing wheels is chopped into a stream of pulses as the timing 
wheel rotates. A pair of photomultiplier tubes converts these light 
pulses to electrical pulses which are then counted electronically. 

Thus the speed of the disk is measured. Fiber optic light guides are 
used for the transport of the light beams throughout this system. 

The sliding velocities are calculated as the difference of the measured 
velocities of the disks. 

Filtered oil is supplied to the exit side of the contact allowing 
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the lubricant one revolution with the disk before entering the con- 
junction. The lubricant is pumped from a 5 gallon supply tank, with 
a thermostatically controlled electric immersion heater and circulator, 
at rates up to 1 gallon per minute. The filter has a paper filter 
element which removes particles down to 1 micron. 

The surface temperature of the disk as it enters the contact zone 
is monitored by an iron-constantan thermocouple trailing on the moving 
surface. An ice bath reference junction is used with the thermocouple. 
Crook [ 44 ] has demonstrated that this method gives accurate results; 
and Johnson and Cameron T 11 ] have found this method agrees closely 

with the temperatures measured by a thermocouple embedded in the surface 
of the disk. 

The electrical output signals of the strain gauges, photomultiplier 
tubes and the thermocouple are continually monitored by a scanning 
digital voltmeter. 

A surface trace of a disk, shown in Figure 5.6, indicates that 
the disks were manufactured with a maximum peak-to-va lley roughness 
of 4 micro-inches. 

5.2 The Lubricants 

The experimental program consisted of the gathering of extensive 
friction data for the two experimental fluids, Mobil XRM 109 F4 and 
Mobil XRM 177 F4. Mobil XRM 109 F4 is a synthesized paraffinic hydro- 
carbon base fluid. Mobil XRM 177 F4 is comprised of Mobil XRM 109 F4 
formulated to improve its anti-fatigue properties. Table 5.1 is the 
physical property data, kindly supplied by the Mobil Research and 
Development Corporation, determined on Mobil XRM 109 F4. The proper- 
ties of Mobil XRM 177 F4 are expected to be the same within experi- 
mental error. 
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Table 5.1 


PHYSICAL PROPERTIES OF MOBIL XRM 109 F4 



Kinematic Viscosity, cs @ 400°F 6.0 

@ 210°F 40.4 

@ 100°F 447 

@ 0°F 37,000 

Total Acid No. 0.0 

Flash Point, °F 520 

Fire Point, °F 595 

Pour Point, °F -60 

Density @ 100°F 0.8389 

@ 200°F 0.8082 

@ 300°F 0.7777 

@ 400°F 0.7428 

Specific Heat @ 300°F 0.635 

@ 400°F 0.692 

Autogeneous Ignition Temp., F 760 

Surface Tension 30.9 
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5.3 Test Procedure 


The strain gauges, photomultiplier tubes, voltage supplies, 
and digital voltmeter must warm up and reach a stable temperature be- 
fore any calibrations are performed. The lubricant supply is also 
heated to the desired temperature during this time. After the warm-up, 
the strain gauge bridges measuring the normal load and frictional 
torque are calibrated. The air cylinder gauges are calibrated to zero 
load, while the torquemeter gauges are calibrated against a shunt re- 
sistance simulating a known torque. 

The oil supply is then turned on and the disk machine may be 
started at minimum load with the field resistances of the two 40 hp 
machines at equal settings. The load and rolling speed are then in- 
creased to the desired values and the bearing torque is measured at 
pure rolling conditions. 

The sliding speed is now varied, while maintaining a constant 
mean rolling speed, to obtain the data for a friction versus sliding 
speed curve. It is easiest to keep the surface temperature within a 
5 degree C range by making some high slip torque measurements first 
and then returning to the low and middle slip values. 

5.4 Results 

A typical set of experimental results is shown in Figure 5.7. 

The friction coefficient, defined as the tractive force divided by the 
applied normal load, is plotted against the sliding speed. The maxi- 
mum Hertzian pressure, rolling speed and lubricant inlet temperature 
remain constant. The friction coefficient rises from zero at pure 
rolling to a maximum value and then decreases with any further in- 
crease in sliding speed. 
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Figure 5.7. Typical experimental results. 



Complete results of the experimental study are presented with 
discussion in the next chapter. 
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CHAPTER VI 


DISCUSSION OF RESULTS 

A new e lastohydrodynamic friction analysis has been developed 
in Chapters II, III and IV. The separate effects of both shear rate 
and time have been included. Shear viscoelasticity results in a 
non-Newtonian relation between the shear stress and shear rate, while 
compressiona 1 viscoelasticity results in a time -dependent viscosity 
function. A numerical solution of the momentum and energy equations, 
with pressure, temperature and time-dependent parameters, is achieved. 

Traction measurements have been made on two synthesized hydro- 
carbon fluids under eiastohydrodynamic conditions. The experimental 
apparatus and procedure have been described in Chapter V. 

This chapter discusses the results of these analytical and experi- 
mental programs. A good correlation of the friction coefficients 
determined by analysis and experiment is shown. 

6.1 Values of the Friction Coefficient Determined by Experiment 

The tractive force transmitted by a thin lubricant film under 
eiastohydrodynamic conditions has been measured for a wide range of 
loads and sliding speeds at high rolling speeds. Specifically, the 
loads ranged from 115,000 psi maximum Hertzian stress to 250,000 psi; 
the sliding speeds varied from zero to over 60 in/sec ; the high rolling 
speeds were 500 and 1000 in/sec; and the oil entrance temperatures 
were 175 °F and 220 °F. The friction coefficient, or. traction coef- 
ficient, calculated from this data for two synthetic paraffinic fluids, 
show variations similar to those found for other lubricating oils by 
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Johnson and Cameron [ 11 ] , Crook [ 4 ] and Bell, and Kannel and 

Allen [ 6 J • 

The friction coefficient rises to a maximum value with increasing 
sliding speed and then decreases with any further increase in the 
sliding speed. The coefficient is also found to increase with increasing 
pressure and to decrease with increasing. rolling speed and temperature. 
Any parameter variation that results in an increase in the friction 
coefficient also results in the maximum friction occurring at a lower 
sliding speed. Examples of this behavior for both experimental fluids 
are shown in Figures 6.1 through 6.14 where the friction coefficient 
is plotted as a function of sliding speed for fixed values of maximum 
Hertzian pressure, rolling speed and oil inlet temperature. 

This behavior may be explained in terms of the liquid model that 
has been developed. At low values of sliding and therefore, low 
shear rate, there is no appreciable temperature gradient across the 
lubricant film. The shear stress increases with shear rate according 
to the effective viscosity predicted by the compressions 1 viscoelastic 
model developed in Chapter III. At slightly higher sliding speeds , 
the temperature rise in the fluid film is significant and cannot be 
neglected. For the range of conditions under study, the analysis 
predicts a rise of film temperature of 15 °F to 20 °F at the sliding 
speed corresponding to the maximum friction coefficient. For even 
higher sliding speeds, it is hypothesized that the mechanism of flow 
changes and is dominated by a pressure and temperature-dependent 
limiting shear modulus.' The temperature at the center of the lubri- 
cant film at the highest sliding speeds is calculated to be 100 °F 
to 150 °F higher than the surface temperature of the disks. 
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Figure 6.2. The effect of load on the friction coefficient. 



EFFECT OF PRESSURE 

U-500 

■"T=i75 



Figure 6.3. The effect of load on. the friction coefficient 
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Figure 6.4. The effect of load on the friction coefficient. 
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Figure 6.5. The effect of lubricant inlet temperature on the friction coefficient. 
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Figure 6.6. The effect of lubricant inlet temperature on the friction coefficient. 
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Figure 6.8. The effect of rolling speed on the friction coefficient 
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Figure 6.9. The effect of rolling speed on the friction coefficient. 
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Figure 6.10. The effect of rolling speed on the friction coefficient. 
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Figure 6.12. The effect of rolling speed on the friction coefficient. 









Figure 6.15. Simplified limiting shear modulus-pressure relationship. 
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Although the actual limiting shear modulus probably follows the dashed 
line, the linear portion at higher pressure levels may be given by 
the relationship shown. Equation (6.2) therefore yields 

u, a m - — (6.3) 

P 

and the friction coefficient also increases with pressure at. high 
sliding speeds. 

An increase in the inlet temperature of the lubricant results 
in lower values of both the viscosity and limiting shear modulus. 

The entire friction coefficient versus sliding speed curve is thus 
lowered according to equations (6.1) and (6.2). This is seen in Figures 
6.5 through 6.7. 

As the rolling speed increases, the film thickness increases 
according to equation (4.3). At low sliding speeds, this higher film 
thickness will reduce the shear rate, which in turn reduces the shear 
stress according to the viscoelastic fluid model. At higher values 
of sliding speed, an increased film thickness results in a higher film 
temperature as indicated by equation (4.18). This will then lower the 
limiting shear modulus. Either of these results, the lower shear stress 
or the lower limiting shear modulus, causes lower friction coefficients 
as seen in Figures 6.8 through 6.14 for the two experimental lubricants. 

6.2 Correlation of Values of the Friction Coefficient Determined 

J^L Experiment and Analysis 

The values of the friction coefficient determined by experiment 
for the Mobil XEM 109 F4 synthetic paraffinic base fluid are compared 
with those predicted by the new analysis in Figures 6.16 through 6.23. 
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Figures 6.24 through 6.31 show a similar comparison for the Mobil 
XRM 177 F4 paraffinic fluid with anti-fatigue additives. A further 
correlation is shown in Figures 6.32 through 6.37 for some of the 
experimental data of Johnson and Cameron f 11 J . 

As in the previous section, the values of the friction coefficient 
are plotted as a function of the sliding speed U^-U^, for fixed values 
of maximum Hertzian pressure P, rolling speed U and oil entrance tem- 
perature T. The units of P, U and T are psi, in/sec and degrees F, 
respectively, except for the Shell Turbo 33 correlations where the oil 
entrance temperature is given in degrees C. Values determined by 
experiment are shown as data points on the curves, while values pre- 
dicted by the analysis are shown by smooth curves. 

The analysis predicts friction coefficients that show the same 
variations as observed experimentally. The friction coefficients rise 
to a maximum value and then decrease with increasing sliding speed; 
they increase with increasing pressure and decrease with increasing 
rolling speed and oil temperature. 

Good correlation is found between the experimental data for Mobil 
XRM 109 F4 and Shell Turbo 33 and the values predicted by the analysis 
using the straight exponential viscosity function adopted by Cheng 
[ 5 ] and the hyperbolic liquid model with c = .25. This corresponds 

to the Barlow, Erginsav and Lamb [ 23 ] liquid model with a limiting 

shear stress. In most cases, the friction coefficients agree within 
10%, with a few extreme cases differing by less than 25%. 

It is necessary to use the hyperbolic model with c = .20 to 
obtain the same correlation for the Mobil XRM 177 F4 lubricant. This 
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Figure 6.18. Correlation of theoretical friction with experimental data. 
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Figure 6.19. Correlation of theoretical friction with experimental data. 
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Figure 6.20. Correlation ol theoretical friction with experimental data. 
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Figure 6.26, Correlation of theoretical friction with experimental data. 
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Figure 6.27. Correlation of theoretical friction with experimental data. 
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Figure 6.28. Correlation of theoretical friction with experimental data 
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Figure 6.30. Correlation of theoretica 1 friction with experimental data. 
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Figure 6.33. Correlation of theoretical friction with experimental data. 




136 


Figure 6.34. Correlation of theoretical friction with experimental data. 
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Figure 6.35. Correlation of theoretical friction with experimental data. 
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Figure 6.36. Correlation of theoretical friction with experimental data. 




lubricant consists of the Mobil 109 F4 as a base with an anti-fatigue 
polymer additive. This additive may change the limiting shear modulus 
function which corresponds to a change in the hyperbolic model constant 
c. A more likely possibility, however, is that the additive increases 
the film thickness. If this is the case, the shear rates would be 
lower in the low sliding speed region and the film temperatures would 
be higher in the high sliding speed region. As previously discussed 
in section 6.1, this would lower the entire friction coefficient curve. 
Thus, if the additive does cause an increase in the film thickness, 
the hyperbolic model with c = .25 might hold true for this lubricant 
also. 

Figure 6.38 shows a comparison of friction versus sliding speed 
curves analytically determined using the viscosity relationships 
adopted by Cheng (equation 4.27) and Chu and Cameron (equation 4.26). 
The curves are extremely close at low sliding speeds but begin to 
diverge at higher sliding speeds as the temperature rise in the lubri- 
cant film becomes larger. The divergence is due to the higher tempera- 
ture dependence of the Chu and Cameron formulation. It is of little 
consequence which formulation is used at low sliding speeds since the 
compressions 1 viscoelastic effects dampen the effects of small changes 
in equilibrium viscosity. The Cheng formulation gives a slightly 
better correlation with all experimental data. Unitl extremely high 
pressure viscosity data is available for lubricants, there will be no . 
other means of choosing among empirical pressure-temperature-viscosity 
functions. 
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Figure 6.38. Comparison of friction curves based upon different viscosity functions. 



6.3 Fluid Property Profiles 


In addition to calculating the friction coefficient, the numerical 
solution of the momentum and energy equations also determines the 
fluid property profiles in the lubricant film. The profiles confirm 
the qualitative estimates of Plint [ 13 ] . 

As the lubricant enters the contact zone, the temperature, and 
therefore the viscosity and the limiting shear modulus are constant 
across the film. At small sliding speeds and low pressures, the 
temperature across the film remains constant and equal to the disk 
surface temperature. The velocity profile is linear and the other 
property profiles are easily predicted. 

Under more severe conditions such as higher sliding speeds and 
pressures, the thermal effects dominate the profiles. The tempera- 
ture profile becomes parabolic, and at the most severe conditions, 
almost triangular. The central plane temperature is 100 °F to 150 °F 
higher than the disk surface temperature. This results in a sharp 
S-shaped velocity profile with an enormous velocity gradient at the 
central plane of the lubricant film. The viscosity, and usually more 
important under these conditions, the limiting shear modulus have 
minimum values on the central plane. Thus, even though the material 
properties and the fluid flow are continuous, the conditions are close 
to those that would occur in a fluid undergoing a discontinuous shear 
failure on the plane of minimum limiting shear stress. 

Two examples of the profiles at the center of the contact zone, 
those of temperature, velocity, viscosity and limiting shear modulus, 
are shown in Figures 6.39 and 6.40 for the analysis of Mobil XRM 109 F4 
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at 200,000 psi maximum Hertzian pressure, 500 in/sec rolling speed and 
175 °F oil inlet temperature. Figure 6.39 is for a 2 in/sec sliding 
velocity where the velocity profile is no longer linear. Figure 6.40 
is for a 50 in/sec sliding velocity where the shear rates at the center 
plane are extremely high. 

6.4 Effect of Convective Heat Transfer 

The simplified energy balance given by equation (4.18) was derived 
by assuming the heat transported by convection was negligible as 
compared with the heat conducted in the two disks. A ratio of con- 
vection to conduction is estimated in equation (4.17) as 

pcUh 2 
2 bk 

2 

Convection has its largest effect for a maximum value of (Uh /b) . 

This corresponds to the condition of maximum rolling speed and minimum 
load. A program was written to include the effects of convection. 
Figure 6.41, the friction coefficient-sliding speed curves, include 
and neglect the convective heat transfer. 

As expected, the convection will carry some heat from the contact 
zone and the lubricant will be slightly cooler. For example, under 
the conditions for maximum convection, the mid-film temperature at the 
center of the contact zone is determined by the analysis to be 5 °F 
cooler; that is, 210 °F compared with 215 °F, at the point of maximum 
friction coefficient. This only affects the friction coefficient at 
higher sliding speeds where the temperature gradient in the film be- 
comes significant. 

The program including the convective effects takes four times 
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Figure 6.41. The effect of convective heat transfer on the friction coefficient. 



the tine required by the simpler program which considers conduction 
only. The results, differing by 5-10% at a maximum, do not warrant 
this expenditure. 

6,5 Effect of Compressional Viscoelasticity 

The overall effect of compressional viscoelasticity on the 
friction coefficient is seen in Figure 6.42. The most prominent 
feature is the shifting of the maximum value to higher sliding speeds. 
This is the same effect a longer shear relaxation time would have on 
the curve. At higher sliding speeds, the flow is dominated by the 
limiting shear modulus. Therefore, to study the _ effects of compres- 
sional viscoelasticity, attention is focused on the region of low 
sliding speeds. 

The values of traction coefficient for very low sliding speeds 
have been calculated are are shown in Figures 6.43, 6.44 and 6.45 
plotted as a function of the ratio of sliding speed to rolling speed, 

Ef = U^-l^/U, for fixed values of peak pressure and U. To simplify 
comparisons with experimental data, the calculations have been made 
for the conditions used by Johnson and Cameron [ 11 1 ; that is, 

hard steel disks of 3 in diameter, a lubricant of viscosity 84 cP at 
atmospheric pressure and 30 °C, and maximum Hertzian pressures of 87, 
110, 147, 176 and 224 x 10 3 psi. 

The variation of traction coefficient with rolling speed may also 
be presented in terms of an effective viscosity. The effective vis- 
cosity r) is defined as that constant viscosity which, for a contact 
area of width 2b and uniform thickness 2h, would give rise to the 
measured tractional force. Then the effective viscosity may be 
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Figure 6.42. The effect of congressional viscoelasticity on the friction coefficient. 
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Figure 6.45. Traction coefficient versus slide/roll ratio, at peak Hertzian pressures of 
200,000 psi and 224,000 psi, for different rolling speeds. 




calculated from the traction coefficiently the expression 

- _ nWfa 

11 " < u rV b 

Alternatively, 



(6.4) 


(6.5) 


where g,/E is the initial slope of the traction curve when plotted as a . 
function of the slide/roll ratio. Values of -T) calculated in this way, 
in the limit of zero sliding speed, are plotted as a function of rolling 
speed in Figure 6.46. 

The calculated values of traction coefficient show a small dependence 
on rolling speed when plotted as a function of the ratio of sliding 
speed to rolling speed. Johnson and Cameron [ 11 ] report that the 

traction coefficient was experimentally found to depend only on the 
slide/roll ratio and to be independent of the rolling speed. In the 
present analysis, this behavior is found only at the lower pressures, 
and then only over a limited range of rolling speeds (Figure 6.43), 
although the values of traction coefficient are similar to those meaured 
by Johnson and Cameron. 

The curvature of the lines in Figures 6.43, 6.44 and 6.45 reflects 
the departure from Newtonian behavior of the lubricant with increasing 
shear rate. At low sliding speeds, the heat generated due to shearing ; 
in the lubricant is negligible. No temperature gradient exists within 
the lubricant film, and the temperature throughout the contact zone 
remains equal to the disk temperature. The decreasing slope of the 
traction coefficient curves with increasing sliding speed is thus a 
consequence of the decrease in the apparent viscosity with increasing 
shear rate. The effect is most marked at the higher pressures and lower 
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Peak Pressure (psixlO 3 ) 



rolling speeds when the small values of film thickness result in' higher 
values of shear rate. The variation of traction coefficient with sliding 
speed is shown in more detail in Figure 6.47 for a peak pressure of 
176,000 psi. The behavior at the other measures conforms to the same 
general pattern. In Figure 6.47 Newtonian behavior, a viscosity which 
is independent of shear rate, is shown by a straight line of unity 
slope. 

It may be calculated from Figure 6.47 that at low rolling speeds 
it is experimentally impossible to obtain Newtonian conditions, as the 
low sliding speeds required -- less than 0.01 in/sec -- are well below 
the experimental range. This fact has important consequences when 
attempts are made to evaluate the effective viscosity from experimental 
data at lower rolling speed. 

Two features of Figure 6.46, showing the variation of effective 
viscosity with rolling speed under isothermal and Newtonian conditions, 
merit special attention. The first of these is the great similarity 
in the shapes of the curves at rolling speeds above 50 in/sec, at all 
but the lowest pressure. This type of behavior has been observed 
experimentally, as seen, for example, in Figure 6.48 taken from Crook 
[4 ]. The second feature is the extremely rapid fall in the effective 

viscosity at low rolling speeds for pressures above 110,000 psi. Re- 
liable measurement of the traction force is difficult at rolling speeds 
below 10 in/sec and high values of peak pressure, as the small film 
thickness becomes comparable with the dimensional irregularities of 
the disk surfaces, and full e lastohydrodynamic conditions no longer 
exist. Extrapolation of results obtained at higher rolling speeds is 
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Figure 6.47. Traction coefficient versus sliding speed, at a peak Hertzian pressure of 
176,000 psi, for different rolling speeds. 




Figure 6.48. 


The effective viscosity (Jj. „) as a function of rolling speed, (a) 30 °C; (b) 45 °C. 

— o — deduced from experimental results; , calculated from visco-elastic hypothesis. 

Load = 7-4 x 10 7 dyncm _1 . . 

Curve from Crook f 4 ] . 
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therefore necessary if a value of effective viscosity at zero rolling 
speed is desired. Furthermore, the present analysis predicts that 
shear rate effects will be significant even at the lowest sliding speeds 
which can be reached experimentally. The measured values of the traction 
force will therefore be less than the values which would be obtained 
under Newtonian conditions. This effect, taken in conjunction with 
the rapid change in value of the effective viscosity at low rolling 
speeds, makes the extrapolation of experimental data to zero rolling 
speed subject to extremely large errors, the magnitude of the error 
increasing as the peak pressure is increased. 

It is suggested therefore, that the observation of Johnson and 
Cameron, [ 11 , Figure 15] whereby the same reduction in effective 

viscosity with rolling speed was observed at all pressures, is a 
consequence of the errors inherent in such an extrapolation. If this 
is so, it follows that the sharp change in the rate of increase of 
viscosity with pressure at pressures above 110,000 psi shown in Figure 
6.49, from' Johnson and Cameron's paper [ 11 ] , is also a consequence 

of the errors in extrapolation, and is not a true property of the 
lubricant . 

To explore this possibility in detail, hypothetical values of 
effective viscosity at zero rolling speed have been obtained by extra- 
polation of the curves of Figure 6. 46 ■, ignoring .the calculated values 
at rolling speeds below 50 in/sec. The shear rate dependence of the 
lubricant viscosity is included by using values of effective viscosity 
calculated at a sliding speed of 0.02 in/sec instead of at the limit 
of zero sliding speed. The values so obtained are shown in Figure 6.50, 
plotted as a function of rolling speed over the range 50 to 500 in/sec. 
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Figure 6.49. Measured variation of effective viscosity versus' peak 
Hertzian pressure, from Johnson and Cameron in'] 



Effective Viscosity, rj (Ibf sec in" 2 ) 



Peak Pressure (psixiO 3 ) 






The curves are then extrapolated below 50 in/sec to obtain a hypothetical 
value of effective viscosity at zero rolling speed. These values are 
plotted in Figure 6.51 to a base of peak pressure (dashed line), and 
show a large deviation from the true values of -effective viscosity 
(solid line) at pressures above 120,000 psi. The experimentally de- 
termined values of effective viscosity shown in Figure 6.49 are also 
plotted in Figure 6.51. The close agreement between the hypothetical 
viscosity curve and the experimental values strongly supports the 
contention that the change in the slope of the viscosity-pressure 
curve at high pressures is an artifact arising from the difficulties 
inherent in the extrapolation procedure, and is not a physical pro- 
perty of the lubricant. 

It has been found that by plotting the data of Figure 6.46 on a 

logarithmic, instead of a linear, scale of rolling speed, the separate 

curves for the different pressures can be combined into a single 

normalized curve. The effective viscosity values are normalized with 

respect to the value at the limit of zero rolling speed riy_Q, and the 

* 

rolling speed values are normalized with respect to U , the rolling . 

speed at which the effective viscosity is equal to 0.5r|y_Q. The re- 

suiting curve of log(r)/r)y_Q) versus log(U/U ) is shown in Figure 6.52. 

* 

The variation of U with the peak pressure in the contact ,is shown in 
Figure 6.53. These two graphs provide a quick and simple method of 
determining the variation of the effective viscosity with rolling speed 
for a given value of maximum pressure. 

In this study of the role of compressiona 1 viscoelasticity in a 
rolling contact system, it has been necessary to simplify the analysis 
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Figure 6.51. Effective viscosity versus peak Hertzian pressure 

; values at zero rolling speed, in the limit 

of zero sliding speed 

; values obtained by extrapolation to zero 

rolling speed, from Figure 6.50 
O ; measured values, from Johnson and Cameron [ \\ 
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Figure 6.52. Normalized plot of effective viscosity versus rolling 
speed . 
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as developed in Chapter III. Among these simplifications, the most 
important are the use of a viscoelastic model with only a single 
retardation time, and the assumption of a symmetrical viscosity dis- 
tribution over the contact zone. It has also been necessary to esti- 
mate the viscosity and the bulk modulus as discussed in detail 
in section 3.3. The simplifications could be eliminated in a more 
detailed analysis. Such improvements are of little va lue , however , 
unless they are matched by improved information about the physical 
properties of the lubricant under the extreme conditions found in 
bearings and other heavily loaded contacts. 

6.6 Comparison of Thermal Theories 

A comparison of several thermal analyses is shown in Figure 6.54 
These include the Johnson-Crook analysis and the author’s numerical 
analysis for the Maxwell- limiting shear stress model and the hyper- 
bolic liquid model (c = .25) corresponding to the B. E. L. -limiting 
shear stress model. Friction coefficients. were also calculated 
neglecting the effects of compressiona 1 viscoelasticity. It is ap- 
parent that these effects must be considered to predict an accurate 
value of the maximum friction coefficient at the correct value of 
sliding speed. 

6.7 Summary and Conclusions 

1. The results of this study demonstrate that values of friction 
coefficient calculated according to the hyperbolic liquid model 
(c = .25) have a good correlation with those determined by experi- 
ment for the two lubricants , .Mobil XRM 109 F4 and Shell Turbo 33. A 
similar correlation is obtained using c = .20 for Mobil XRM 177 F4. 
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Figure 6.54. Comparison of the experimental friction curve for Shell Turbo 33 at P = 176,000 psi, 
U = 260 in/sec and T = 30 °C and those calculated by several theoretical analyses. 


This change in the value of c might possibly be explained by an in- 
creased film thickness due to the polymer additives in this fluid. 

This hypothesis is described in section 6.2. 

2. The effects of shear rate and time are separated and explained by 
the two phenomena of shear viscoelasticity and compressiona 1 visco- 
elasticity, respectively. 

3. A unified description of the non-Newtonian shear rate dependence 
of the viscosity is presented as a new hyperbolic liquid model. With 
this model, the transition from the non-linear region to the shear 
modulus dominated region is shown to be a smooth one. In the high- 
slip region, where the friction is dominated by the shear modulus, 
the variation of friction with load is very sensitive to the pressure 
dependence of the shear modulus. 

4. The friction coefficient rises to a maximum value with increasing 
sliding speed and then decreases with any further increase in the 
sliding speed. The coefficient is also found to increase with in- 
creasing load and to decrease with increasing rolling speed and tem- 
perature . 

5. The effects of compressiona 1 viscoelasticity are developed in terms 
of a simple model for the volume creep of a liquid following the ap- 
plication of a pressure step. This model is used to determine the de- 
pendence on rolling speed of the friction coefficient between highly 
loaded rolling contacts. Curves are presented which show the vari- 
ation with rolling speed of the effective viscosity of the lubri- 
cant in the contact zone under isothermal conditions. Both the shape 
of the curves and the values of effective viscosity are consistent 
with the results of experimental measurements. The shape of the curves 


166 



in this region is found to be nearly independent of the peak pressure 
in the contact. 

6. At very low values of rolling speed, in a region which is experi- 
mentally inaccessible, the analysis predicts a very rapid variation 
of effective viscosity with rolling speed. It is shown that, as a 
consequence, the extrapolation of experimental data to zero rolling 
speed can result in extremely large errors in the estimated values of 
effective viscosity. 

The results of this study suggest future work that will increase 
the understanding of friction in elastohydrodynamic lubrication. The 
most urgently needed research is in the field of fluid rheology. The 
viscosity and density of lubricants at high pressures would be extremely 
helpful, and shear and compressional relaxation experiments must be 
performed to measure the fluid moduli at high pressures. This work 
is needed to confirm and expand our understanding of the mechanism 
of flow under EHD conditions. 
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APPENDIX A 

NUMERICAL INTEGRATION OF A 
NON-EQUIDISTANTLY ' TABULATED FUNCTION 


A method of overlapping parabolas is employed to yield a second 
order approximation to the integral of a non-equidistant ly tabulated 
function graphically represented in Figure A.l. > 

The function f(x) may be represented by a second order Taylor’s 


series expansion about x^: 

(x-x ) 

i(x) = f(x ) + — f*(x ) + 

n il n 


(x-x 

— f"(x 1 + 

2! U n' 


Accordingly, 


f = f + h f ’ .+ f " + • 
r»+l n nn 2 n 


(A. 2) 


f = f - h f’ + f" + 

n-1 n n-1 n 2 n 


(A. 3) 


f n +1 

* £(x n + l 

f n-1 

s £< Vl> 
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n 
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X 
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1 


Equations (A.2) and (A. 3) are solved simultaneously to yield 

2 2 2 

f n = h T(h " + h ) f n-l " h A T f n + h (h . + h ) f n+ 1 (A ' 4) 

n-1 n-1 n' n-1 n n n-1 n 


-h (h -h .) h . 

f i _ E f + — 2 — 2— L_ f Ull f 

n h ,(h . + h ) n-1 h n h n h (h . + h ) n + 1 

n-1 n-1 n n-1 n n n-1 n 


(A. 5) 


168 



f(x) 



Figure A.l. Graphical representation of the non-equidistantly tabu- 
lated function f(x). 
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The shaded area of Figure A.l is calculated by integrating f(x) 

-between the limits x and x The function is approximated by a 

n n + 1 

parabolic curve through f ,, f and f This integral is called 

n- 1 n n + I 

T n + 1 

i 1 : 

1 n 


n + i P n+1 

l*n “ J f(x) dx 

x 

n 




dx 


2 3 

h h 

= h f +T S f' 
n n 2 n 


+ TT~ f" 

6 n 


(A. 6) 


The derivatives are evaluated by equations (A. 4) and (A. 5), resulting 
in the following expression: 


T n+ 1 
In 


-h‘ 


6h . (h . + h ) 
n-1 n-1 n 


f i + 
n-1 


h (3h + h ) 

n n-1 n 


6h 


n-1 


h (3h . + 2h ) 

n n-1 n_ 

6(h + h ) 

n-1 n 


n + 1 


(A. 7) 


Similarly, the same integral may be calculated using a parabolic 
approximation for f(x) through the points f^, f^ + and f Equa- 
tion (A. 1) determines 


f = f + (h + h )f ’ + |(h , . + h ) 2 f * 
n + 2 n n + 1 n n * n+1 n n 


(A. 8) 


This second integral, 

pansion about x . , : 
n+1 


is defined by the Taylor's series ex- 
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(A. 9) 


x 

, in+1 =f | £ 4 . 1 + ^ X " X + l'( X " X , , ,1 

l n j |_ n + 1 n+ln+1 * n+1 n + lj 


= h f - — £ 1 + — f * 

n n+1 2 n + 1 6 n+1 


The derivatives f" , , and f* , , are evaluated in exactly the same 
n+1 n + 1 

manner as determined equations (A. 4) and (A. 5). These results are sub- 
stituted into equation (A. 9) yielding 

„ . . h (2h + 3h .. ) h (h + 3h . .) 

n + 1 n n n+1 f n n n + 1 f 

2 n _ 6(h + h . n + 6h . , n + 1 

n n+1) n+1 


+ 6h 1 (h + h ) f n+2 
n+1 n n+1 


The best possible second order approximation of the integral is 
the average of the two values just calculated; this is the method of 
overlapping parabolas. Therefore, the average integral, I^ + is 
defined as 


I n+1 = 1 [ I n+1 +• I n + 1 l 
n ^ [1 n 2 n J 


where and 2^n + ^ are cie ^^ net ^ by equations (A. 7) and (A. 10) , 

respective ly . 

The integral of the function f(x) between the limits x = a and 
x = b is defined by equation (A. 12): 


£(*) dx „ ji; + 1 + 1 


n=a + 1 


In order to keep the integral a function only of values in the 
range of integration, the method of overlapping parabolas is not used 


on the two extreme intervals. 


171 



APPENDIX B 


NUMERICAL ANALYSIS 


The FORTRAN IV coding of the friction analysis is listed in 
this appendix. It is followed by examples. of data cards and the 
optional subroutines. 
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PROGRAM CONTROL! INPUT, TAPE60= INPUT « OUTPUT t PUNCH ) CONTR 1 

♦ ♦♦♦♦♦♦♦♦♦♦♦♦♦■fr********************* CONtR 2 





♦CONtR 

3 

* 

PROGRAM CONTROL 


♦CONTR 

A 

* . 



♦CONTR 

5 

* 



♦CONTR 

6 

* 

CALCULATES THE FRICTION AND VELOCITY AND TEMPERATURE PROFILES 

♦CONTR 

7 

* 

IN AN ELASTOHYDRODYNAM IC LUBRICATED CONTACT. 


♦CONTR 

8 

* 



♦CONTR 

9 

* 

REQUIRED SUBPROGRAMS - 


♦CONTR 

10 

* 

SUBROUTINE PRINTS 


♦CONTR 

1 1 

* 

FUNCTION VISC 

- 

♦CONTR 

12 

* 

FUNCTION ZERO 


♦CONTR 

13 

■x 

FUNCTION PS I 


♦CONTR 

1 A 

* 

SUBROUTINE SECANT 


♦CONTR 

15 

* 

FUNCTION EXP I 


♦CONTR 

16 

* 

SUBROUTINE I NTEG 


♦CONTR 

17 

* 

SUBROUTINE RTM I 


♦CONTR 

18 

X 



♦CONTR 

19 

* X -X 


********* 

♦ CONTR 

20 

* 



♦CONTR 

21 

tt 

PROGRAM SET UP FOR MOBIL DATA NEGLECTING 

CONVECTION 

CONTR 

22 


COMMON C.GC.LLLL 


CONTR 

23 


COMMON /CPSI/ GInF.ETA0«T«Y*H*U2U1 «NP,DUDY, IC* 

TG. OMEGA. CH* AH 

CONTR 

2A 


COMMON /CPSIO/ XK1.XK0 


CONTR 

25 


COMMON /CO/ Q 


CONTR 

26 


COMMON /CZERO/ I PTRANS * TRANS ( 6 I 


CONTR 

27 


COMMON /CPR/ ETA2 


CONTR 

28 


DIMENSION G INF (21 I *ETAO (21 )«T(21)*Y<21> *DUDY( 21 ) *TG( 21 ) ,0MEGA(21 ) CONTR 

29 


DIMENSION Q( 21 ) 


CONTR 

30 


DIMENSION C(2«3«21 ) * GC (21*21 ) 


CONTR 

31 


DIMENSION ET A2 ( 2 1 ) 


CONTR 

32 


DIMENSION DTDY ( 2 1 ) *NEU/T(21 1 .TRACTCF(20> ,SL IP(20> 

CONTR 

33 


EXTERNAL PS I 


CONTR 

34 


REAL logeta*newt 


CONTR 

35 


INTEGER COUNT 


CONTR 

36 

* 



CONTR 

37 

•* 

PHYSICAL CONSTANTS AND DATA 


CONTR 

38 

* 



CONTR 

39 

* 

IRDP MAXIMUM HERTZIAN PRESSURE / 1000. 

- , 

CONTR 

40 

* 

I RDU MEAN ROLLING SPEED 


CONTR 

41 

* 

IRDT LUBRICANT INLET TEMPERATUPE 


CONTR 

42 

* 

PH IT THERMAL REDUCTION FACTOR 


CONTR 

43 

* 

AH, CH HYPERBOLIC LIQUID MODEL PARAMETERS 


CONTR 

44 

9000 

READ 9001 « IRDP* IRDU« I RDT « PH I T « AH , CH 


CONTR 

45 

9001 

FORMAT (4X13* 15. 14 *F6.3* 10XF3. 1 «6XF4.2> . 


CONTR 

46 


IF (EOF (60 ) >9999,9002 

- 

CONTR 

47 

9002 

CONTINUE 


CONTR 

40 


PUNCH 9003 


CONTR 

49 

9003 

FORMAT!/) 


CONTR 

50 


PHZ=!000.*FL0AT( IRDP) 


CONTR 

51 


U=E LOA T ( I RDU) 


CONTR 

52 


TOIL=FLOAT< I RDT ) +460 • 


CONTR 

53 

* 

NP NUMBER OF GRID POINTS ACROSS THE HALF-FILM THICKNESS 

CONTR 

54 


NP=11 


CONTR 

55 


NH = NP - 1 


CONTR 

56 

* 

NIP NUMBER OF PRESSURE STEPS IN HALF CONTACT-LENGTH 

CONTR 

57 


NIP = 3 


CONTR 

58 


173 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


* 


* 

* 

* 


TW = TO I L 
TMAX=950. 

COND THERMAL CONDUCTIVITY OF THE LUBRICANT 

COND= • 1 *778./3600. 

CYLW CONTACT WIDTH OF THE DISKS 

CYLW=0.25 

Rl. R2 RADII OF THE DISKS 
Rl =R2=3.0 

El* E2 ELASTIC MODULUS OF THE DISKS 
El=E2=30.E+6 

POIS1* POIS2 POISSONS RATIO FOR THE DISKS 
PO I S 1 =PO I S2=0 • 3 

ALPHA VISCOSITY PRESSURE COEFFICIENT FOR THE LUBRICANT 

ALPHA= 1 . 04E-4 
BET A =5 • 1E7*ALPHA 
GAMMA = 930 * * ALPHA 

DK THERMAL CONDUCTIVITY OF THE DISKS 

DK = 21 . 7*778. /3600. 

DRHO DENSITY OF THE DISKS 

DRHO= • 283 

DC SPECIFIC HEAT OF THE DISKS 

DC=. 109*778. *12. 

HERSA* HERSB CONSTANTS FOR THE HERSHEL VISCOSITY EQUATION 
HERS A=8 .974 
HERSB=-3 .2 

NGRAPH REQUIRED NUMBER OF GRAPHS FOR EACH TEMPERATURE PROFILE 
NGRAPH=0 

MGRAPH REQUIRED NUMBER OF GRAPHS OF TRACTION COEF . VS SLIP 
MGRAPH= 1 

PRNT = 2HON GIVES ADDITIONAL OUTPUT FOR DEBUGGING PURPOSES 

PRNT = 3H0FF 

INITIALIZATION AND BOUNDARY CONDITIONS 

COUNT = 0 
IC = 1 

DTDY ( 1 )=0»0 

Y< 1 ) =0.0 

TRACT = 0.0 

FLASH=0. 0 

PI =3.1415927 

DO 10 1=1 «NP 

NEWT ( I ) = 1 O'. * < 1 . — Y < I ) )+TW 

IF (I.LT.NP) Y(I+1)=Y(I)+1. /FLOAT ( NH ) 

10 CONTINUE 

R=R 1 *R2/ < R 1 +R2 ) 

E = 2./( < 1 .-POIS1*POIS1 J/El + I 1 . —PO I S2*PO I S2 ) /E2 > 

B=4.*R/E*PHZ 

FLASHK=0 .24 /SORT < P I *PHZ*DK*DRHO*DC *U*R/E ) 

TRANS! 1 ) = R/U/6 . 

TRANSI2) =B/U/6. 

TRANS < 3) =B/U/2. 

CALCULATION OF load 

W=PHZ*PHZ*PI*(2./E)*R 
PRINT 2 * PHZ ♦ W ♦ R * E 

2 FORMAT (* 1 CONTROL PHZ = *E15.8»* W = *E15.8** R = *E15.8* 


CONTR 59 
CONTR 60 
CONTR 61 
CONTR 62 
CONTR 63 
CONTR 64 
CONTR 65 
CONTR 66 
CONTR 67 
CONTR 68 
CONTR 69 
CONTR 70 
CONTR 71 
CONTR 72 
CONTR 73 
CONTR 74 
CONTR 75 
CONTR 76 
CONTR 77 
CONTR 78 
CONTR 79 
CONTR 80 
CONTR 81 
CONTR 82 
CONTR 83 
CONTR 84 
CONTR 85 
CONTR 86 
CONTR 87 
CONTR 88 
CONTR 89 
CONTR 90 
CONTR 91 
CONTR 92 
CONTR 93 
CONTR 94 
CONTR 95 
CONTR 96 
CONTR 97 
CONTR 98 
CONTR 99 
CONTR 100 
CONTR 101 
CONTR 102 
CONTR 1 03 
CONTR 104 
CONTR 105 
CONTR 106 
CONTRl 07 
CONTR 108 
CONTR 1 09 
CONTRl 10 
CONTRl 1 1 
CONTRl 12 
CONTRl 13 
CONTRl 14 
CONTRl 15 
CONTRl 16 
CONTRl 17 
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* 

* 


* 


* 


* 

* 


* 

* 


,* 


* 


* 

* 


* 


* 


1 * E = *E15.8) 

CALCULATION OF HALF-FILM THICKNESS 

ET AENT = 1 0 4 ** ( HERS A+HERSB*ALOG 10(TOIL— 460* > )*1 .45E-7 

H= 1 .6*ALPHA**0.6*< ETA£NT*U) **0 4 7*E**0.03*R**0.43/W**0. 13 

H= 1 .2*H 

H=PHIT*H 

H=0.5*H - 

PRINT 1*H«ALPHA,ETAENT*U«CH4AH 
1 FORmAT(*OH = *E12.5«* ALPHA = *E12*5.* ETAENT = *E12.5. 

1* U = *F 6* 0 ♦ * CH = *F5«3»* AH = *F7.3) 

SLIP LOOP 

NSL I P= 15 

DATA (SLIP(IU)«IU=l«20)/.5«l.*2.«3.«4««5»t6.«8.«10»*15.»20. 
1 « i 50« « 6C • 4 5*0 • 0/ 

DO 6000 I U = 1 4 NSL I P 

DATA (TRACTCFI IP) 4 IP=1 4 20) /20#0 4 0 / 

FLASH=FLASHK*TRACT*SLIP< IU) 

PRINT 6 4 FLASH 

6 FORMAT ( * CONTROL FLASH = *E15.8) 

U2U 1=0 .5# SL I P ( I U ) 

PRINT 84 IU.SLIPI IU) 

8 FORMAT <*1 CONTROL IU = *134* SLIP = *E15.8> 

HERTZIAN PRESSURE LOOP 

TRACT = 0 « 0 
DO 5999 IP= 1 .NIP 
IPTRANS= IP 

XB=(2.*FL0AT< IP)-1. ) /2 • /FLOAT ( N I P ) 

P=PHZ*SQRT ( XB* ( 2.-XB ) ) 

SOLVE' MOMENTUM EQUATION 

I TCOUNT — 0 

4 IF ( PRNT • EQ # 2H0N ) PRINT 44 4 COUNT 
44 FORMAT <*OCOUNT = *13) 

COUNT = 0 
DO 11 1= 1 4NP 

IF ( ITCOUNT.EQ.O.OR. ITCOUNT.GT. 10 ) T(I)=NEWT(I) 

IF (IT COUNT 4 GT • O « AND • I T COUNT 4 LE 4 1 0 ) T < I ) =0 .5* ( T ( I > +NEWT ( I ) ) 
IF ( I TCOUNT 4 GT . 1 00 ) GO TO 6003 
IF (T(I).GTiTMAXI T(I)=TMAX 
TW=TOIL+FLASH 

GlNFt I )= 1 .2*P/(2.52+.01333*(T( I )-492. > )-l .45E4 
IF ( G I NF ( I ) 4 LT 4 1 4 ) GINF(I) = 1. 

CALCULATION OF STEADY-STATE VISCOSITY 
ETEXP=ALPHA*P+ ( BET A + GAMMA*P > * ( 680 • — T ( I ) )/680.XT( I ) 

ET A2 ( I ) = 4 62*EXP (ETEXP)*1 4 45E-5 
CALCULATION OF TIME-DEPENDENT VISCOSITY 
ET AO ( I ) = V I SC ( P 4 ET A2 ( 1)4 IVCODE) 

11 CONTINUE 

DUMMY=PSI0(0.0) 

CALL SECANT (XKO 4 XK1 4 T AU 4 PS I 4 XK 4 4 00 1 4 500 4 I TERR > 

IF ( I TERR • EO 4 1 ) STOP 
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CONTRl 
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33 
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43 
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50 
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• IF ( PRNT • EQ • 2HON ) PRINT 41.TAU 
41 FORMAT!*OTAU = *E15.8) 

* 

* SOLVE ENERGY EQUATION 

* 

DO 5000 I = 1 ,NP 
Q( I ) =-TAU*DUDY( I )/C0ND*H 
5000 CONTINUE 

IF (PRNT.EQ.2H0N) CALL PR 1 
DO 5010 1=2* NP 

CALL INTEG1 (0. *Y( I ) « Y«Q*NP*DTDYC I > *C,GC.LLLL, IERR) 

5010 IF ( IERR.NE.O) PRINT 5012* IERR 
5012 FORMAT < * INTEGRATING O* IERR = *13) 

DO 5020 1=2* NP 

CALL INTEG1 CO. ,Y( I ) « Y . DTDY * NP « NEWT ( I) *C*GC*LLLL* IERR) 
5020 IF! IERR. NE.O) PRINT 5022*!ERR 
5022 FORMAT ( * INTEGRATING DTDY* IERR = *13) 

XK7=TW-NEWT(NP) 

I T COUNT = I TCOUNT + 1 
NEWT ( 1 ) = 0 • O 
DO 5030 1=1, NP 

NE WT ( I ) = NEWT ( I ) +XK7 

5030 ; I FI ABSI 1 « — I T ( I ) /NEWT! I))).LT. 0.001) COUNT = COUNT+ 1 
IFICOUNT.NE.NP) GO TO 4 
IF I IVCODE.NE.O) PRINT 5034 
.5034 F0RmAT(*0 MINIMUM VISCOSITY REDUCTION*) 

CALL PR 1 
PRINT 99 

PRINT 5040. I (NEWT! I ) * I ) * 1 = 1 , NP > 

5040 FORMAT!* NEW TEMP = *F7.2** I = *14) 

* 

* CALCULATION OF TRACTION 

* 

TRACT = TRACT + 2. *B/N I P*TAU 

5999 CONTINUE 

TRACTCF ( I U ) = TRACT /W 
PRINT 6002* TRACTCF I IU ) 

6002 FORMAT (*OCONTROL TRACTCF! IU) = *E15.8> 

* 

* PLOT TEMPERATURE PROFILE 

* 

. IF (nGRAPH. EQ.O) GO TO 6000 
DO 5052 NGR=1 .NGRAPH 
, PRINT 5050, SLIP! IU) 

5050 FORMAT I # 1 (U2 - Ul) = *E15.8> 

PRINT 99 

CALL STPL Till, NEWT , Y , NP * 1 H* , 1 , 1 HY > 

PRINT 5051 

5051 FORMAT! 1H0.92X, 1 1 HTEMPERATURE ) 

5052 CONTINUE 

6000 CONTINUE 

6003 PRINT 98 

PRINT 6001 . ! ITRACTCF! IU) *SL IP! IU) ) * IU=1 .NSLIP) 

6001 FORMAT!* CONTROL TRACT = *E15.8,* U2-U1 = *E15.8>. 

PUNCH 6005* ( (SLIP! IU) .TRACTCF! IU) ). IU=1 *NSLIP) 

6005 FORMAT ( 2F I 0 » 5 ) 

* 

* PLOT TRACTION COEFFICIENT 


CONTR176 
CONTR 1 77 
CONTR 1 78 
CONTR 179 
CONTR 180 
CONTR 181 
CONTR 182 
CONTR 183 
CONTR 184 
CONTR 185 
CONTR 186 
CONTR 187 
CONTR 188 
CONTR 189 
CONTR 190 
CONTR 191 
CONTR 192 
CONTR 193 
CONTR 194 
CONTR 195 
CONTR 1 96 
CONTR 197 
CONTR 198 
CONTR 199 
C.0NTR200 
CONTR201 
C ONTR202 
C0NTR203 
CONTR204 
CONTR205 
C0NTR206 
CONTR207 
C0NTR208 
C0NTR209 
CONTR210 
CONTR21 1 
CONTR2 1 2 
C0NTR213 
C0NTR214 
C0NTR2 1 5 
C0NTR2 1 6 
C0NTR2 1 7 
C0NTR2 1 8 
C0NTR2 19 
CONTR220 
C0NTR221 
CONTR222 
CONTR223 
CONTR224 
C0NTR225 
C0NTR226 
CONTR227 
CONTR228 
CONTR229 
C0NTR230 
C0NTR23 1 
C0NTR232 
CONTR233 


* 


C0NTR234 
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IF CmGRAPH.EO.O) GO TO 6013 
DO 6012 MGR= 1 tMGRAPH 
PRINT 98 

CALL STPLT1 < 1 ,SL I P , TRACTCF ,NSL 1 P « 1 H* . 14 , 14HTRACT ION COEF.) 
PRINT 60 1 1 

60 11 FORMAT ( 1 HO, 100X*7HU2 - Ul) 

6012 CONTINUE 

6013 CONTINUE 

98 FORMAT! *1*1 

99 FORMAT(*0*> 

GO TO 9000 

9999 CONTINUE 
STOP 
END 


C0NTR23S 
CONTR236 
C0NTR237 
C0NTR238 
CONTR239 
C0NTR240 
CONTR24 1 
CONTR2A2 
CONTR2A3 
CONTR244 
CONTR245 
CONTR2A6 
C0NTR2A7 
C0NTR2 A 8 
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SUBROUTINE PRINTS PRINT 

COMMON /CPSI/ G I NF (ETAUiTtYiHt U2U 1 * NP • DUD Y . IC 1 TG 1 OMEGA 1 CH 1 AH PRINT 

COMMON /CQ/ Q PRINT 

COMMON /CP R/ ETA2 PRINT 

DIMENSION GINFI21 ) «ETA0(21 ) ,T (21 > « YI21) .DUDY (21) «TG(21 ) * OMEGA ( 2 1 ) PRINT 
DIMENSION Q { 2 1 ) PRINT 

DIMENSION ET A2 ( 2 1 ) PRINT 

ENTRY PR 1 PRINT 

PRINT 1 PRINT 

1 FORMAT <-*0 I Y < I > T(I) G I NF t I ) TG ( 1 ) ETAO(I) PRINT 

1 ETA2 ( I ) OMEGA! I) OUDY(i) 0 ( 1 )*) PRINT 

PRINT 99 PRINT 

PRINT 3 . ( C I » Y< I ) «T( I 1 «GINF( I ) tTG( I ) «£TAO< I ) * ET A2 ( I ) » OMEGA C I ) .DUDY (PR I NT 
1 I ) «Q< I ) ) . 1 = 1 .NP) PRINT 

3 FORMAT** *I2.F7.2.F9.2.7<2XE1 1 .4) ) PRINT 

99 FORMAT(*0*> PRINT 

RETURN PRINT 

END PRINT 


1 

2 

3 

4 

5 

6 

7 

8 
9 

1 0 
1 1 
1 2 
1 3 

14 

15 

16 
1 7 
18 
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FUNCTION VISC!P.ETA2«C0DE> 




VISC 

1 

* 

* * 

********************** 

* * 

* * 

* * 

* * * * * VISC 

2 



* 




♦vise 

3 

* 


FUNCTION VISC 




♦vise 

4 

* 






♦vise 

5 



* 




*V I sc 

6 

* 


CALCULATES THE TRANSIENT VALUE OF THE LUBRICANT 

VISCOSITY. *VISC 

7 



* 




*vl sc 

8 

* 


ARGUMENTS - 




*visc 

9 

* 


P PRESSURE 




♦vise 

1 0 

* 


ETA2 EQUILIBRIUM VALUE OF VISCOSITY 




*VISC 

1 1 

* 


CODE ERROR PARAMETER 




*VISC 

1 2 

* 






*VISC 

13 

* 


REQUIRED SUBPROGRAMS - 




*visc 

1 4 

* 


FUNCTION ZERO 




*visc 

15 

* 


FUNCTION EXP I 




*V I sc 

1 6 

* 






*VISC 

1 7 

* 


COMMON STORAGE - 




*VISC 

18 

* 


THE VARIABLE TAUP MUST BE IN LABELED COMMON 

CV1SC. 


*visc 

19 

* 






*V!SC 

20 

* 


ERROR INDICATIONS - 




*VISC 

21 

* 


CODE = 0 INDICATES NO ERROR. 




*VI sc 

22 

* 


CODE = 1 INDICATES THE TRUE TRANSIENT VISCOSITY IS 

OUT 

OF THE *VISC 

23 

* 


RANGE OF THE PROGRAM AND THE MAXIMUM 1 

POSSIBLE 

*visc 

24 

* 


VISCOSITY REDUCTION WAS ASSUMED. 




*visc 

25 

* 






*visc 

26 

* 

* * 

*** ** ***************** 

* * 

* * 

* * 

* * * * * VISC 

27 



* 




*visc 

28 



COMMON /CVISC/ TAUP 




visc 

29 



EXTERNAL ZERO 




visc 

30 



REAL KF 




VISC 

3 1 



INTEGER R TM I ERR . CODE 




VISC 

32 



PRNT=3H0FF 




VISC 

33 



CODE=0 




VISC 

34 



KF=3.5E5+9.*P 




VISC 

35 



TAUP = 50 . #ET A2/KF 




VISC 

36 



IF (PRNT.EO.2HCN) PRINT 81. KF « T AUP » P* ETA2 




VISC 

37 


ai 

FORMAT!* VISC KF = *E13.6«* TAUP = *E13 

• 6 « * 

P 

= *E13.6* VISC 

38 



1* ETA2 = *E 13.6) 




VISC 

39 

* 






V I SC 

40 

* 


TEST OF RANGE 




vise 

4 1 

* 

■, 





vise 

42 



IF (ZERO! -.0 1 >*ZERO< -20. ) .GT.O. ) GO TO 60 




VISC 

43 

* 






VISC 

44 

* 


ROUGH BOUNDING OF ZERO 




VISC 

45 

* 






VISC 

46 



IF ( ZERO ( — 5 » ) *ZERO (—20. ) .GT • 0 » ) GO TO 30 




VISC 

47 

* 


ROOT IS BETWEEN -20. AND -5« 




VISC 

48 



XL I =-20. 




VISC 

49 



XRI =-5. 




VISC 

50 



GO TO 50 




VISC 

51 

* 


ROOT IS BETWEEN -5. AND -.01 




VISC 

52 


30 

I COUNT =0 




VISC 

53 



XL I = — 5 . 




VISC 

54 



XR I =— 3 . 




VISC 

55 



ZEROXP I = ZERO ( XR I ) 




VISC 

56 



ZEROXL I =ZERO ( XL I ) 




VISC 

57 


40 

IF ( ZEROXR I * ZEROXL I » LE . 0 » ) GO TO 50 




VISC 

58 
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IF (XRI .GT.-.Ol ) GO TO 60 

VISC 

59 



XL I = XR 1 

VISC 

60 



XR I =XR I /3 « 

VISC 

61 



ZEROXL I =ZEROXR I 

VISC 

62 



ZEROXR I = ZERO ( XR 1 ) , 

vise 

63 



I COUNT = I COUNT + 1 

VISC 

64 



IF C ICOUNT.GT. 100) GO TO 92 

vise 

65 



GO TO 40 

VISC 

66 

-* 



. VISC 

67 

*• 


PRECISE DETERMINATION OF ZERO 

VISC 

68 

* 



VISC 

69 


50 

CALL RTM l ( S I » ZEROS I • ZERO « XL I * XR I « • 05. 1 00 * RTM I ERR ) 

VISC 

70 



IF (RTMIERR.NE.O) GO TO 90 

vise 

71 



VISC=ETA2*EXP!SI ) 

VISC 

72 



RETURN 

VISC 

73 

-* 



VISC 

74 

* 


OUT OF RANGE 

VISC 

75 

* 



VISC 

76 


60 

IF (ZEROC-20. ) .lT.O. ) GO TO 70 

. VISC 

77 

* 



VISC 

78 

* 


NEGLECT VISCOSITY REDUCTION OF LESS THAN 1. PER CENT 

vise 

79 

-* 



VISC 

80 



VISC=ETA2 

VISC 

81 



RETURN 

vise 

82 

* 



vise 

83 

-* 


MAXIMUM VISCOSITY REDUCTION 

VISC 

84 

* 



visc 

85 


70 

VISC=ETA2*EXPf-20. ) 

VISC 

86 



CODE= 1 

VISC 

87 



RETURN 

VISC 

88 

* 



vise 

89 

* 


ERROR MESSAGES 

VISC 

90 

* 



VISC 

91 


90 

PRINT 91 * RTM I ERR 

VISC 

92 


91 

FORMAT ( * vise ERROR in RTMI ERR = *12) 

VISC 

93 



STOP 

VISC 

94 


92 

PRINT 93 

VISC 

95 


93 

FORMAT!* VISC MAX NO. OF ITERIONS EXCEEDED IN PRERTM I 

PROC.*) VISC 

96 



STOP 

VISC 

97 



END 

VISC 

98 
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FUNCTION ZERO (SI) 



ZERO 

I 

* * * 


* 

* * * 

* ZERO 

2 

* 




♦ ZERO 

3 

-* 

function zero 



♦ ZERO 

A 

-* 




♦ ZERO 

5 

* 




♦ ZERO 

6 

* 

THIS EXTERNALLY SUPPLIED FUNCTION IS NEEDED BY SUBROUTINE 

RTMI 

♦ ZERO 

7 

* } 

CALLED FOR IN FUNCTION VISC. IT MYST BE PRESENT WHEN A 

TIME- 

♦ZERO 

8 

* 

dependent viscosity function is used. 



♦ ZERO 

9 

* 




♦ ZERO 

10 

* * * 


* 

* * * 

* ZERO 

1 1 


COMMON /CZER O/ IPTRANS.TRANSI6) 



ZERO 

12 


COMMON /CVISC/ TAUP 



ZERO 

1 3 


INTEGER EXPI ERR 



ZERO 

14 


PRNT = 3H0FF 



ZERO 

IS 


IF (SI.GE.-20.) GO TO 10 



ZERO 

16 


SI--20. 



ZERO 

17 


PRINT 20 



ZERO 

18 

20 

FORMAT (*ZERO SI HAS BEEN READJUSTED TO -20.*) 



ZERO 

19 

' 10 

CONTINUE 



ZERO 

20 


CALL EXPI (EISI .Si *EXPIERR) 



zero 

2 1 


IF (EXPIERR.NE.O) GO TO 90 



ZERO 

22 


ZERO=TRANS( I P TRANS I +T AUP* < E I S I —0 »05*EXP t S I > 1 



ZERO 

23 


IF (PRNT.EQ.2H0N) PRINT 8 I « S I . E I S I * TRANS I I PTRANS ) • ZERO 



ZERO 

24 

81 

FORMAT ( * ZERO SI =*E13.6.* EISI = *E13.6»* TRANS 

= 

*E 13. 

6 .ZERO 

25 


1* ZERO = *E 13.6) 



ZERO 

26 


RETURN 



zero 

27 

*. . 




ZERO 

28 

* 

ERROR MESSAGE 



ZERO 

29 

* 




ZERO 

30 

90 

PRINT 91 .EXP I ERR 



ZERO 

31 

91 

FORMAT(* ZERO ERROR IN EXPI ERR = *13) 



ZERO 

32 


STOP 



ZERO 

33 


END 



ZERO 

34 
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FUNCTION PSItTAU) 


PS I 

1 

* * 

■a-****-***********-***-**** 

**•**#**** 

* PSI 

2 




*PSI 

3 


FUNCTION PS I 


♦PSI 

4 




*PS I 

5 




*PS I 

6 


THIS FUNCTION SUBROUTINE SUPPLIES THE RHEOLOGICAL MODEL 

*PS I 

7 


FOR THE LUBRICANT. IN THIS CASE. WE ARE USING 

THE 

*PS I 

8 


hyperbolic MODEL. 


*PS I 

9 




*ps r 

10 

* * 

•****•*-****■**■*-«■•***-***•**** 

**•**■*##** 

* PS I 

1 1 


COMMON C.GC.LLLL 


PS I 

12 


COMMON /CPSI/ GInF.ETAO.T. Y.H.U2U1 .NP.DUDY. IC. 

TG, OMEGA. CH. AH 

PS I 

13 


COMMON /CPSIO/ XK1.XK0 


PSI 

14 


DIMENSION G I NF (21 1 *ETA0(21 ) «T <21 > «Y(21 1 .DUDYI21 > «TG(21 ) . OMEGA (21 ) PS I 

15 


DIMENSION C(2*3*21). GC (21.21) 


PSI 

16 


REAL intg 


PSI 

17 


PRNT=3H0FF 


PS I 

18 

100 

CONTINUE 


PSI 

19 


IF (PRNT.EQ.2H0N) PRINT 200.TAU 


PS I 

20 

200 

FORMAT(*OTAU = *E 15.8) 


PSI 

21 


DO 210 1=1, NP 


PSI 

22 


TG< I ) = T AU/G I NF ( I ) 


ps r 

23 


IF(TG( I ) .LT.CH) GO TO 205 


PS I 

24 


TAU=0.99999999*CH*GlNF ( I ) 


PS I 

25 


GO TO 100 


PS I 

26 

205 

OMEGA ( I ) =CH*TG( I )*( AH*CH*TG( 11 + (TGI I >-CH>*(TG( I )-CH) )/(CH~TG( I ) >**PSI 

27 

13 


PS I 

28 


DUD Y ( I ) = G I NF < I l/ETAOf I >*OMEGA{ I ) *H 


PSI 

29 

210 

CONTINUE 


PS I 

30 


GO TO (220.230) . IC 


PSI 

31 

220 

CALL INTEG ( O . , 1 , , Y • DUD Y . NP « I NTG . C . GC , LLLL . I ERR ) 

PSI 

32 


IC = 2 


PS I 

33 


GO TO 240 


PSI 

34 

230 

CALL I NTEG2 ( 0 . , 1 ♦ .Y.DUDY.NP, I NTG . C , GC » LLLL ♦ I ERR > 

PS I 

35 

240 

IF (IERR.NE.0) PRINT 241.IERR.TAU 


PSI 

36 

241 

FORMAT ( * I ERR = *14,* AT TAU = *E15.8) 


PS I 

37 


PS I = ALOG 1 C ( I NTG/U2U 1 ) 


PS I 

38 


IF (RRNT.EQ.2H0N) PRINT 250, INTG, PS I 


PS I 

39 

250 

FORMAT ( * INTG = *E15.8,* PSI = *E15.8) 


PS I 

40 


TDUDY=0. 0. 


PSI 

41 


DO 310 1=1, NP 


PSI 

42 

310 

TDUDY = TDUDY +DUD Y ( I ) 


PSI 

43 


AD J=U2U1 *FLOAT ( NP ) /TDUDY 


PSI 

44 


IF ( AOJ.GT.0.99) GO TO 390 


PS I 

45 


DO 320 1=1. NP 


PS I 

46 

320 

DUD Y ( I ) = AD J*DUD Y ( I ) 


PS I 

47 

390 

CONTINUE 


PSI 

48 


RETURN 


PS I 

49 


ENTRY PS 10 


PS I 

50 


XK 1 =0 » 20*CH*G I NF ( 1 ) 


PSI 

51 


XKO=0 . 5*XK 1 


PSI 

52 


PS I =0.0 


PSI 

53 


RETURN 


PS I 

54 


END 


PSI 

55 
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SUBROUTINE SEC ANT ( XO . X 1 * XF I NAL « FUNC tX i CONV « MAX 1 T * I ERR > SECNT 1 

♦ ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦* SECNT 2 






♦SECNT 

3 

* 


SUBROUTINE SECANT 


♦SECNT 

4 

* 




♦SECNT 

5 





♦SECNT 

6 

* 


SECANT METHOD SOLUTION OF FUNC(X) = 0, 

♦SECNT 

7 

-* 

-*** 

ALLOWS POSITIVE X ONLY 


♦SECNT 

8 

* 




♦SECNT 

9 

* 


ARGUMENTS - 


♦SECNT 

10 

* 


XO* XI TWO INITIAL GUESSES OF 

ROOT 

♦SECNT 

1 1 

* 


XF I NAL FINAL ESTIMATE OF ROOT 


♦SECNT 

12 

* 


FUNC EXTERNALLY SUPPLIED FUNCTION FUNC 

♦SECNT 

13 

* 


X PARAMETER OF FUNCTION 

FUNC 

♦SECNT 

14 

* 


CONV TEST FOR CONVERGENCE 


♦SECNT 

15 

-8- 


MAX IT MAXIMUM NUMBER OF ITERATIONS TO FIND SOLUTION 

♦SECNT 

16 

* 


I ERR ERROR PARAMETER 


♦SECNT 

17 

* 




♦SECNT 

18 

* 


REQUIRED SUBPROGRAMS - 


♦SECNT 

19 

* 


FUNCTION FUNC ( X ) 


♦SECNT 

20 

* 




♦SECNT 

21 

* 


COMMON STORAGE - NONE 


♦SECNT 

22 

* 




♦SECNT 

23 

* 


ERROR INDICATIONS - 


♦SECNT 

24 

* 


I ERR = 0 INDICATES NO ERROR. 


♦SECNT 

25 

* 


I ERR = 1 INDICATES THE MAX. NO. 

OF ITERATIONS were exceeded 

♦SECNT 

26 

* 




♦SECNT 

27 

* 


EDWARD G. TRACHMAN 

M.E. DEPT. 492-5640 

♦SECNT 

28 

* 




♦SECNT 

29 

* 

* * 

tt*-*-*****-*-**-*--#--**-* 


* SECNT 

30 



PRNT = 3H0FF 


SECNT 

31 



IT = 0 


SECNT 

32 



I ERR=0 


SECNT 

33 



XMIN=-1 . E99 


SECNT 

34 



XMAX= 1.E99 


SECNT 

35 



FXMIN=-2.E99 


SECNT 

36 



FXMAX= 2.E99 


SECNT 

37 



FO-FUNC ( XO ) 


SECNT 

38 



IF (FO.LT.O.) XMIN=X0 


SECNT 

39 



IF (FO.LT.O.) FXMIN=FO 


SECNT 

40 



IF (FO.GT.O.) XMAX=XO 


SECNT 

41 



IF (FO.GT.O.) FXMAX=FO 


SECNT 

42 



FI =FUNC(X1 ) 


SECNT 

43 


1 

IF (FI .LT.O. .AND. XI .GT.XMIN) GO 

TO 50 

SECNT 

44 


51 

IF (Fl.GT.O.. AND .X1.LT. XMAX ) GO 

TO 52 

SECNT 

45 


53 

IF (PRNT.EO.2HON) PRINT 2 * XO * X 1 


SECNT 

46 


2 

FORMAT** SECANT XO = *E 1 5 • 8 * * 

XI = *E 15.8) 

SECNT 

47 



IF (ABS(Fl).LT. CONV ) GO TO 10 


SECNT 

48 



IT = IT + 1 


SECNT 

49 



IF ( I T • GT • MAX IT) GO TO 99 


SECNT 

50 



SLOPE = ( F 1 — F 0 ) / ( X 1 -XO ) 


SECNT 

51 



CLOSE=DIM(ABS<FG) . ABS(F1 ) > 


SECNT 

52 



IF (PRNT.EQ.2H0N) PRINT 80*X0*X1 *F0«F1 « SLOPE 

SECNT 

53 


80 

FORMAT! *OSECANT XO = *E15.8** 

XI = *EI5.8,* FO = *E 15.8.* 1 

FI SECNT 

54 


1= *E15.8,* SLOPE = *E15.8) 


SECNT 

55 



IF (CLOSE. GT.. 0001 ) 5.6 


SECNT 

56 

* 


USE XI 


SECNT 

57 


cc 

DELX=-F1 /SLOPE 


SECNT 

58 
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X0 = X1 
F0 = F1 
GO TO 7 
* USE XO 

6 DEL X=—FO /SLOPE 

7 XI =XO+DELX 

*** ALLOWS POSITIVE X ONLY 
IF (Xl.LE.O.) X 1 = 1 • E— S 
IF (Xl.GT.XMIN) GO TO 8 
X1=XMIN 
F1=FXMIN 
GO TO 1 

8 IF (Xl.LT.XMAX) GO TO 9 
X1=XMAX 

F1=FXMAX 
.GO TO 1 
- 9 FI =FUNC ( XI) 

GL TO 1 
to XFInAL=X1 

IF (PRNT.EQ.2H0N) PRINT 11. XF INAL .Ft* IT 
11 FORMAT ( * G SECANT XF I nAL = *E15.B«* F(XFINAL) = *E15.8i* IT 

1*14.) 

RETURN 
50 XM I N=X 1 
FXM I N=F 1 
GO TO 51 
52 XM AX -XI 
FXMAX=F1 
GO TO 53 
99 PRINT 100 

100 FORMAT!* SECANT MAXIMUM NUMBER OF I TERAT I ONS EXCEEDED*) 

IERP= 1 
RETURN 
END 


SECNT 59 
SECNT 60 
SECNT 61 
SECNT 62 
SECNT 63 
SECNT 64 
SECNT 65 
SECNT 66 
SECNT 67 
SECNT 68 
SECNT 69 
SECNT 70 
SECNT 71 
SECNT 72 
SECNT 73 
SECNT 74 
SECNT 75 
SECNT 76 
SECNT 77 
SECNT 78 
= SECNT 79 
SECNT 80 
SECNT 81 
SECNT 82 
SECNT 83 
SECNT 84 
SECNT 85 
SECNT 86 
SECNT 87 
SECN T 88 
SECNT 89 
SECNT 90 
SECNT 91 
SECNT 92 
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SUBROUTINE EXPI (REStXt IERR) 




EXPI 

1 

* * 

-a--**-***-**-**-*****-*-**-*-****** 

* * 

* 

* * * * * 

* EXPI 

2 

* 





♦ EXPI 

3 

* 

SUBROUTINE EXPI 




♦ EXPI 

4 

* 





♦EXPI 

5 

* 





♦EXPI 

6 

* 

COMPUTES THE EXPONENTIAL INTEGRAL FOR NEGATIVE 

ARGUMENTS. 

♦ EXPI 

7 

* 

IN THE RANGE -20 TO ZERO . 




♦ EXPI 

8 

* 





♦ EXPI 

9 

* 

FOR X EQUAL TO 0 THE RESULT VALUE IS SET TO 1.E75. 



♦EXPI 

10 

* 

FOR X LESS THAN -20 OR GREATER THAN ZERO THE CALCULATION IS 

♦EXPI 

1 1 

* 

BYPASSED AND THE ARGUMENT REMAINS UNCHANGED. 




♦ EXPI 

12 

* 





♦EXPI 

13 

* 

THE EXPONENTIAL INTEGRAL IS DEFINED AS THE 




♦ EXPI 

14 

* 

RES INTEGRAL ( ExP ( -T ) /T . SUMMED OVER T FROM X TO 

INFINITY) • 

♦ EXPI 

15 

* 





♦ EXPI 

16 

* 

A POLYNOMIAL APPROXIMATION IS USED FOR ARGUMENTS IN 

THE 

♦EXPI 

17 

* 

RANGE -5 TO ZERO. 




♦EXPI 

18 

* 

REF. LUKE AND WIMP, -JACOBI POLYNOMIAL EXPANSIONS 

OF 

A 

♦ EXPI 

19 

* 

GENERALIZED HYPERGEOMETRIC FUNCTION OVER A SEMI 

-INFINITE RANGE- 

. ♦EXPI 

20 

* 

mathematical tables And other aids to COMPUTATION. 



♦ EXPI 

21 

* 

VOL. 17. 1963. ISSUE 84, PP« 395-404. 




♦ EXPI 

22 

* 





♦ EXPI 

23 

'* 

AN EXPONENTIAL APPROXIMATION IS USED FOR ARGUMENTS 

IN 

THE 

♦EXPI 

24 

* 

RANGE -20 TO -5. 




♦EXPI 

25 

.* 





♦ EXPI 

26 

* 

arguments - 




♦ EXPI 

27 

* 

RES RESULT VALUE. 




♦EXPI 

28 

* 

X ARGUMENT OF EXPONENTIAL INTEGRAL. 




♦EXPI 

29 

* 

IERR RESULTANT ERROR PARAMETER. 




♦ EXPI 

30 

* 





♦ EXPI 

31 

* 

REQUIRED SUBPROGRAMS - NONE 




♦EXPI 

32 

* 





♦ EXPI 

33 

* 

COMMON STORAGE - NONE 




♦EXPI 

34 

* 





♦ EXPI 

35 

* 

ERROR INDICATIONS - 




♦EXPI 

36 

* 

IERR a 0 INDICATES NO ERROR. 




♦ EXPI 

37 

* 

IERR a 1 INDICATES X IS LESS THAN -20. 




♦EXPI 

38 

* 

IERR a 2 INDICATES X IS POSITIVE. 




♦EXPI 

39 

* 





♦ EXPI 

40 

* 

EDWARD G. TRACHMAN M.E. 

DEPT. 

492-5640 

♦EXPI 

41 

* * 

■fr******************.***.*** 

* * 

* 

***** 

♦ EXPI 

42 

* 





EXPI 

43 

* 

TEST OF RANGE 




EXPI 

44 

* 





EXPI 

45 


IERR=0 




EXPI 

46 


IF (X.GT.O. ) GO TO 60 




EXPI 

47 


IF (X.GT.-5.) GO TO 35 




EXPI 

48 


IF (X.LT.-20.) GO TO 55 




EXPI 

49 

* 





EXPI 

50 

* 

ARGUMENT IS BETWEEN -20 AND -5. 




EXPI 

51 

* 





EXPI 

52 

* 

C=-l .094 14E-3/3.**-5 




EXPI 

53 


C=— 2.65876020E— 01 




EXPI 

54 


RES=C*3.**X 




EXPI 

55 


RETURN 




EXPI 

56 

* 





EXPI 

57 

* 

ARGUMENT IS BETWEEN -5 AND ZERO. 




EXPI 

58 
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35 


* 


* 


* 

* 


* 

* 


* 


* 

* 

* 


40 


x=-x 

I F ( X 1 40 4 50 *40 

PES = -ALOG < ABS ( X ))- ( (((((<<((<((• 1 03 17602E- 1 1 *X“. 1 5798675E- 1 0 > *X+ 
1 . 1 6826592E-9 ) *X- • 2191 5699E-8 ) *X+ • 27635830E-7 ) *X- .3072622 1 E-6 1 *X+ 
2. 3099604 OE-5 ) *X-. 28337590E-4 ) *X+« 2314 83P2E-3)*X-. 00 16666906 >*X+ 
3.01 04 16662 )*X-. 055555520 )*X+. 25 ) *X- 1 . 0 1 *X- . 5772 1 566 


RES=-RES 

x=-x 

RETURN 


ARGUMENT IS EQUAL TO 7ER0. 

50 RES= 1 • E75 
RES=~RES 

x=-x 

RETURN 

ARGUMENT IS LESS THAN -20. 

55 I ERR= 1 
RETURN 

ARGUMENT IS POSITIVE. 

60 I ERR=2 
RETURN 
END 


EXP I 

59 

EXP I 

60 

EXP I 

61 

EXP I 

62 

EXP I 

63 

EXP I 

64 

EXP I 

65 

EXP I 

66 

EXP I 

67 

EXP I 

68 

EXP I 

69 

EXP I 

70 

EXP I 

71 

EXP I 

72 

EXP I 

73 

EXP I 

74 

EXP I 

75 

EXPI 

76 

EXP I 

77 

EXPI 

78 

EXPI 

79 

EXPI 

80 

EXPI 

81 

EXPI 

82 

EXPI 

83 

EXPI 

84 

EXPI 

85 

EXPI 

86 
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* 

* 

* 

* 

-H- 

* 

-* 

* 

* L 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

■* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* * 


* 

* 

* 

* 

* 

* 


SUBROUTINE I NT EG (A«B*X*P« NP « VALUE »Ci GC • L * I ERR ) 

♦ ♦♦♦♦*♦♦. ******** •- 5 r '*^».&*********.**** 

SUBROUTINE INTEG 


INTEGRATES THE NON EQUIDISTANTUY TABULATED FUNCTION FIX!!)) 
BETWEEN THE LIMITS A AND B« WHERE A OR B MUST EQUAL F<X<1)>» 

A MODIFIED METHOD OF OVERLAPPING PARABOLAS IS EMPLOYED. 

ENTRY POINTS - 

INTEG FIRST TIME SUBROUTINE IS CALLED AND 

WHEN NEW ABSCISSAS ARE USED. 

INTEG1 WHEN NEW LIMITS OF INTEGRATION ARE USED. 

INTEG2 WHEN USING THE SAME ABSCISSAS AND 

LIMITS OF INTEGRATION AS THE LAST CALL. 

ARGUMENTS - 

A LOWER LIMIT OF INTEGRATION. 

B UPPER LIMIT OF INTEGRATION. 

X ARRAY OF ARGUMENT VALUES. MUST BE MONOTON I CAlLY 

INCREASING AND MUST BE DIMENSIONED NP. 

F ARRAY OF FUNCTION VALUES. MUST BE DIMENSIONED NP. 

NP NUMBER OF POINTS. NP MUST BE GREATER THAN 3. 

VALUE RESULTANT VALUE OF THE INTEGRATION. 

C. GC WEIGHTING FUNCTIONS PASSED TO THE MAIN PROGRAM 

FOR STORAGE. 

L LIMITS OF INTEGRATION PASSED TO THE MAIN PROGRAM 

FOR STORAGE. 

I ERR RESULTANT ERROR PARAMETER. 

REQUIRED SUBPROGRAMS - NONE 

COMMON STORAGE - c 

THE WEIGHTING FUNCTIONS C AND GC ARE STORED IN THE MAIN PROGRAM 
AND REQUIRE THE FOLLOWING DIMENSION STATEMENT WHERE D. GE .NP. 
DIMENSION C ( 2 » 3. D ) iGC(DtO) 

ERROR INDICATIONS - 


INTEG 1 
* INTEG 2 
♦INTEG 3 
♦INTEG 4 
♦INTEG 5 
♦INTEG 6 
♦INTEG 7 
♦INTEG 8 
♦INTEG 9 
♦INTEG 10 
♦INTEG 11 
♦INTEG 12 
♦INTEG 13 
♦INTEG 14 
♦INTEG 15 
♦INTEG 16 
♦INTEG 17 
♦INTEG 18 
♦INTEG 19 
♦INTEG 20 
♦INTEG 21 
♦INTEG 22 
♦INTEG 23 
♦INTEG 24 
♦INTEG 25 
♦INTEG 26 
♦INTEG 27 
♦INTEG 28 
♦INTEG 29 
♦INTEG 30 
♦INTEG 31 
♦INTEG 32 
♦INTEG 33 
♦INTEG 34 
♦INTEG 35 
♦INTEG 36 
♦INTEG 37 
♦INTEG 38 
♦INTEG 39 
♦INTEG 40 


I ERR = 0 INDICATES NO ERROR. 

I ERR = 1 INDICATES NP IS LESS THAN 4. 

I ERR a 2 INDICATES THE LIMITS OF INTEGRATION ARE NOT At NODES 
OR ARE OUT OF THE RANGE OF THE TABLE. 

EDWARD G. TRACHMAN M.E. DEPT. 492-5640 

♦ ♦♦♦**♦♦♦*■»♦♦***************** *** 
DIMENSION X ( NP) ,F ( NP > .C( 2. 3«NP) ,GC ( NP.NP ) 

DIMENSION H< 100) 

NP MUST BE GREATER THAN 3 

IF (NP.LE.31 GO TO 96 

CALCULATION OF INTERVALS OF X 

NH=NP— 1 

187 ’ 


♦INTEG 41 
♦INTEG 42 
♦INTEG 43 
♦INTEG 44 
♦INTEG 45 
♦INTEG 46 
♦INTEG 47 
* INTEG 48 
INTEG 49 
INTEG 50 
INTEG 51 
INTEG 52 
INTEG 53 
INTEG 54 
INTEG 55 
INTEG 56 
INTEG 57 
INTEG 58 



DO 10 1 = 1. NH 
10 H( I )=X( 1+1 ) —X ( I ) 

DO 20 I = 1 * NH 

IF (I »EQ • 1 ) GO TO 15 

* . 

* DEFINE COEFFICIENTS OF FIRST PARABOLA 

* 

C < 1 « 1 « I ) = -<H { I ) ) ** 3 / ( 6 • *H ( I -1 )*<H< 1-1 >+H( I ) ) ) 

C( 1 .2. I ) =H( I )*.<3.*H< 1-1 )+H< I > )/(6.*H( 1-1 ) ) 

C< 1 .3. I > =H( I >*(3.*H< i-1 )+2.*H( I ) )/(6.*<H< 1-1 >+H( I )) ) 
15 CONTINUE 

IF ( I .EQ.NH) GO TO 20 

* DEFINE COEFFICIENTS OF SECOND PARABOLA 

•* 

C ( 2 * 1 « I > =H< I >*( 2,*H( I )+3.*H ( I + 1 ) ) /< 6.* ( H< I )+H( 1+1 ) ) ) 
C(2 .2 * I ) =H( I >*(H< I )+3.*HC 1 + 1) 1 + 1 1 ) 

C(2«3»I)=-<H(I> )**3/<6.*H( 1+1 )*(H( I>+H( 1+1 ) ) } 

20 CONTINUE 


* DEFINE GROUPED COEFFICIENTS 

* 

DO 61 L=liNP 
DO 6 1 I = 1 « NP 
61 GC ( I » L ) =0.0 

GC ( 1 ♦ 1 ) =C <2. 1 . 1 > 

GC (.2 « 1) =C (2*2. 1 ) 

GC ( 3 . 1 }=C(2.3, 1 ) . 

NPM2=NP— 2 

DO 65 L = 2 • NPM2 

LP2=L+2 

DO 65 1=1. LP2 

C A = 0 • 0 

CB = 0 • 0 

IF ( I-L+2 .GT. 0. AND. I -L+2 .LT.4 ) CA=C(1.I — L+2 .L > 
I F ( I-L+l . GT • 0 • AND « I -L+ 1 .LT.4) CB=C (2. I-L+l »L) 

65 GC f I .L> =GC( I .L-l >+0.5*<CA+CS) 

NPM3=NP-3 

DO 66 1=1. NPM3 

66 GC ( I .NP- 1 ) =GC ( I »NP-2 ) 

GC (NP-2.NP- 1 ) =GC <NP-2 .NP-2 ) +C ( 1 . 1 . NP-1 ) 

GC ( NP- 1 . NP- 1 ) =GC ( NP- 1 « NP— 2 ) +C ( 1 » 2 . NP- 1 ) 

GC ( NP.NP- 1 ) =GC ( NP.NP-2 ) +C ( 1 .3 .NP- 1 ) 

-* 

ENTRY INTEG1 

-* 

* SETTING LIMITS OF INTFGRATION 

* 

IF (B-A) 40.92.30 

* 

* B IS GREATER THAN A 

* 

30 ALIM = A 
BLIM = B 
S I GN = I . 0 
GO TO 50 

* A is GREATER THAN B 

* 


INTEG 59 
INTEG 60 
INTEG 61 
INTEG 62 
. INTEG 63 
INTEG 64 
INTEG 65 
INTEG 66 
INTEG 67 
INTEG 68 
INTEG 69 
INTEG 70 
INTEG 71 
INTEG 72 
INTEG 73 
INTEG 74 
INTEG 75 
INTEG 76 
INTEG 77 
INTEG 78 
INTEG 79 
INTEG 80 
INTEG 81 
INTEG 82 
INTEG 83 
INTEG 84 
INTEG 85 
INTEG 86 
INTEG 87 
INTEG 88 
INTEG 89 
INTEG 90 
INTEG 91 
INTEG 92 
INTEG 93 
INTEG 94 
INTEG 95 
I NTEG 96 
INTEG 97 
INTEG 98 
INTEG 99 
I NTEG 1 00 
I NTEG 101 
I NTEG 102 
I NTEG 103 
I NTEG 104 
I NTEG 105 
I NTEG 106 
I NTEG 107 
I NTEG 108 
I NTEG 1 09 
I NTEG 1 10 
I NTEG 1 1 1 
I NTEG 1 12 
I NTEG 1 13 
INTEGl 14 
INTEG 1 15 
INTEGl 16 
INTEGl 17 
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* 

* 

* 

* 


* 


* 


* 

* 


* 


* 

* 


* 


* 

«■ 

* 

* 


40 ALIM = B 
BLIM a A 
SIGN =-1*0 
50 NH=NP-1 
11=0 

DO 59 I=liNP 

IF ( I I —2 I 55.57,599 

55 XXXX=ALIM-X< I ) 

IF (ABS(XXXX).LT. *0000000001) GO TO 56 
IF (XXXX*GT. 0.0) 57,97 

56 IALIM = I 
11=11+2 

57 IF (II.EQ.l) GO TO 59 
XXXX=BL I M-X ( NP+ 1 - I ) 

IF (ABSIXXXX) .LT.. 0000000001 ) GO TO 58 
IF CXXXX.GT. 0.0) 97.59 

58 I BLIM = NP+l-I 
11 = 11 + 1 

59 CONTINUE 

599 IF (II.NE.3) GO TO 97 - 

IF ( IALIM.NE. 1 ) GO TO 97 
L= IBLIM-1 

ENTRY INTEG2 

CALCULATION OF INTEGRAL OVER SUB INTERVAL 

VALUE = 0.0 
LP2=L+2 

IF CLP2.GT.NP) LP2=NP 
DO 80 1=1 »LP2 

80 VALUE a VALUE +GC (I*L)*F(I) 

CALCULATE THE FINAL VALUE OF THE INTEGRAL 
VALUEaS I GN*VALUE 

SET ERROR PARAMETER FOR NORMAL RETURN 

92 I ERR a 0 
RETURN 

SET ERROR PARAMETER FOR TOO FEW POINTS 

96 I ERR a 1 
RETURN 

SET ERROR PARAMETER FOR A AND/OR B NOT AT NODES 
OR OUT OF RANGE OF THE TABLE 

97 I ERR = 2 
RETURN 
END 


INTEGl 18 
INTEG 1 19 
INTEG120 
INTEG121 
INTEG 122 
I NTEG 123 
INTEG 124 
I NTEG 1 25 
I NTEG 126 
I NTEG 127 
I NTEG 128 
I NTEG 129 
I NTEG 1 30 
I NTEG 131 
I NTEG 132 
I NTEG 133 
I NTEG 134 
I NTEG 135 
I NTEG 136 
I NTEG 137 
I NTEG 138 
INTEG 139 
I NTEG 140 
INTEG141 
I NTEG 142 
INTEG 143 
I NTEG 144 
I NTEG 145 
I NTEG 146 
I NTEG 147 
I NTEG 148 
INTEG 149 
I NTEG 150 
I NTEG 151 
I NTEG 152 
I NTEG 1 53 
I NTEG 154 
I NTEG 155 
I NTEG 156 
INTEG 157 
I NTEG 158 
I NTEG 159 
INTEG 160 
I NTEG 161 
I NTEG 1 62 
I NTEG 163 
I NTEG 164 
I NTEG 165 
I NTEG 166 
I NTEG 167 
I NTEG 1 68 
I NTEG 169 
I NTEG 170 
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TYPICAL DATA CARDS 


0000000001111111111222222222233333333334444444444555555555566666666667777777777b 
123456 78901 x 23456 7890123456 789012 3456 7 890123456 7890123456 789012345678901234567890 


16 

115 

500 

175 

.81 

AH=0 . 0 

CH=0 .25 

PROG (MOBIL) 

vise 

23J71 

16 

154 

500 

175 

.81 

AH=0 . 0 

CH=0. 25 

PROG (MOBIL) 

vise 

23J71 

16 

200 

500 

175 

.81 

AH=0.0 

CH=0 . 25 

PROG (MOBIL) 

vise 

23J71 

16 

250 

500 

220 

.89 

AH=0. 0 

CH=0. 25 

PROG (MOBIL) 

vise 

23J71 

16 

115 

1000 

175 

.58 

AH=0. 0 

CH=0. 25 

PROG (MOBIL) 

vise 

23J71 

16 

154 

1000 

175 

.58 

AH=0. 0 

CH=0. 25 

PROG (MOBIL) 

vise 

23J71 

16 

200 

1000 

175 

.58 

AH=0. 0 

CH=0. 25 

PROG (MOBIL) 

vise 

23J71 

16 

200 

1000 

220 

. 72 

AH=0 . 0 

CH=0 . 25 

PROG (MOBIL) 

vise 

23J71 

16 

115 

500 

175 

.81 

AH=0 . 0 

CH=0 . 20 

PROG (MOBIL) 

vise 

23J71 

16 

154 

500 

175 

.81 

AH=0 . 0 

CH=0 . 20 

PROG (MOBIL) 

vise 

23J71 

16 

200 

500 

175 

.81 

AH=0 . 0 

CH=0. 20 

PROG (MOBIL) 

vise 

23J71 

16 

250 

500 

175 

.81 

AH=0. 0 

CH=0. 20 

PROG (MOBIL) 

vise 

23J71 

16 

250 

500 

220 

.89 

AH=0. 0 

CH=0 . 20 

PROG (MOBIL) 

vise 

23J71 

16 

115 

1000 

175 

.58 

AH=0. 0 

CH=0. 20 

PROG (MOBIL) 

vise 

23J71 

16 

154 

1000 

175 

.58 

AH=0.0 

CH=0. 20 

PROG (MOBIL) 

vise 

23J71 

16 

250 

1000 

220 

.72 

AH=0 . 0 

CH=0. 20 

PROG (MOBIL) 

vise 

23J71 
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PROGRAM CONTROL < I NPUT , TAPE60= I NPUT ♦ OUTPUT ♦ PUNCH ) 

CONV 

1 

* * * 

■a-*#*-*-*-**-*******-*******-********** 

* CONV 

2 

* 


■frCONV 

3 

* 

PROGRAM CONTROL 

*CONV 

4 

* 


*CONV 

5 

* 


♦ CONV 

6 

* 

CALCULATES THE FRICTION AND VELOCITY AND TEMPERATURE PROFILES 

♦ CONV 

7 


IN AN ELASTOHYDRODYNAM IC LUBRICATED CONTACT. THIS VERSION OF 

♦CONV 

8 

* 

THE PROGRAM INCLUDES THE EFFECTS OF CONVECTIVE HEAT TRANSFER. 

♦ CONV 

9 

-* 


♦CONV 

10 

* 

REQUIRED SUBPROGRAMS - 

♦ CONV 

1 1 

* 

SUBROUTINE PRINTS 

♦ CONV 

12 

* 

FUNCTION VISC 

♦CONV 

13 

-* 

function ZERO 

♦CONV 

14 

* 

function PS I 

♦CONV 

15 

* 

subroutine secant 

♦ CONV 

16 

* 

function expi 

♦ CONV 

17 

* 

subroutine integ 

♦CONV 

IB 

-* 

subroutine rtmi 

♦CONV 

19 

* 


♦CONV 

20 

* * * 


* CONV 

21 

* 


♦CONV 

22 

* 

PROGRAM SET UP FOR MOBIL XRM OIL WITH CONVECTION 

CONV 

23 


COMMON C « GC « LLLL 

CONV 

24 


COMMON /CPSI/ G I nF (ETAO iTtYtHi U2U 1 « NP « DUD Y . I C * TG * OMEGA . CH ♦ AH 

CONV 

25 


COMMON /CPSIOX XK 1 * XKO 

CONV 

26 


COMMON /CQ/ Q 

CONV 

27 


COMMON /CZ ERO/ IPTRANS.TRANS(6) 

CONV 

28 


COMMON /CPR / ETA2 

CONV 

29 


DIMENSION GINF ( 21 ) .ETA0I21 ) * T ( 2 1 > . Y ( 2 1 ) . DUDY ( 2 1 > »TGt 21 > .OMEGA t 21 ) CONV 

30 


DIMENSION Q( 21 ) 

CONV 

31 


DIMENSION C ( 2 « 3 « 2 1 ) « GC <21.21 ) 

CONV 

32 


DIMENSION ETA2I21) 

CONV 

33 


DIMENSION DTD Y ( 2 1 ) .NEWT (21 ) .TRACTCFI20) .SLIP! 20) 

CONV 

34 


DIMENSION TJM1 (21 ) 

CONV 

35 


EXTERNAL PS I 

CONV 

36 


REAL LOGETA.NEWT 

CONV 

37 

*■ 


CONV 

38 

* 

PHYSICAL CONSTANTS AND DATA 

CONV 

39 

* 


CONV 

40 

9000 

READ 900 1 . I RDP . IRDU. I ROT « PH I T . AH . CH 

CONV 

41 

9001 

FORMAT (4X1 3. 15. I4.F6.3. 10XF3. 1 .6XF4.2) 

CONV 

42 


IF (EOF (60) ) 9999. 9002 

CONV 

43 

9o02 

CONTINUE 

CONV 

44 


PUNCH 9003 

CONV 

45 

9003 

FORMAT ( / ) 

CONV 

46 


PHZ= 1000 »*FLOAT ( IRDP) 

CONV 

47 


U=FLOAT( I RDU) 

CONV 

48 


T0IL=FL0AT< I RDT ) +460 . 

CONV 

49 

* 

NP NUMBER OF GRID POINTS ACROSS THE HALF-FILM THICKNESS 

CONV 

50 


NP =11 

CONV 

51 


NH = NP - 1 

CONV 

52 

* 

NIP NUMBER OF PRESSURE STEPS IN HALF CONTACT-LENGTH 

CONV 

53 


NIP = 3 

CONV 

54 


NIP2=2*NIP 

CONV 

55 


TW= TO 1 L 

CONV 

56 

* 

COND THERMAL CONDUCTIVITY OF THE LUBRICANT 

CONV 

57 


COND=. 1*778. /360C. 

CONV 

58 

* 

OILRHO DENSITY OF THE LUBRICANT 

CONV 

59 
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OII_RHO=.0325 




CONV 

60 

-* 

OILC SPECIFIC HEAT OF THE LUBRICANT 




CONV 

61 


OILC=.4*778.*12. 




CONV 

62 


CYLW=0.25 




CONV 

63 

* 

CYLW CONTACT WIDTH OF THE DISKS 




CONV 

64 

•H- 

R 1 * R2 RADII OF THE DISKS 




CONV 

65 


R 1 =R2=3.0 ' - 




CONV 

66 

*• 

E I « E2 ELASTIC MODULUS OF THE DISKS 




CONV 

67 


El=E2=30.E+6 




CONV 

68 

* 

POISI* POIS2 POISSONS RATIO FOR THE DISKS 




CONV 

69 


PO I S 1 =PO I S2 = 0 • 3 




CONV 

70 

-* 

ALPHA VISCOSITY PRESSURE COEFFICIENT FOR 

THE 

LUBRICANT 

CONV 

71 


ALPHA= 1 . 04E-4 




CONV 

72 


beta = 5 » i e 7* alpha 




CONV 

73 


GAMMA = 930 » * ALPHA 




CONV 

74 

-fr 

DK THERMAL CONDUCTIVITY OF THE DISKS 




CONV 

75 


DK=21 *7*778. /3600. 




CONV 

76 

* 

DRHO DENSITY OF THE DISKS 




CONV 

77 


DRHO = « 283 




CONV 

78 

* 

DC SPECIFIC HEAT OF THE DISKS 




CONV 

79 


DC= . 1 09*778. *12. 




CONV 

80 

* 

HERSA. HERSB CONSTANTS FOR THE HERSHEL VISCOSITY 

equation 

CONV 

81 


HERS As 8 . 974 




CONV 

82 


HERSBs-3.2 




CONV 

83 

* 

NGRAPH REQUIRED NUMBER OF GRAPHS FOR EACH 

TEMPERATURE 

PROFILE 

CONV 

84 


NGRAPH=0 




CONV 

85 

-* 

MGRAPH REQUIRED NUMBER OF GRAPHS OF TRACTION COEF • VS 

SLIP 

CONV 

86 


MGRAPH = 1 




CONV 

87 

* 

PRNT=2HON GIVES ADDITIONAL OUTPUT FOR DEBUGGING 

PURPOSES 

CONV 

88 


PRNT =3HOFF 




CONV 

89 

* 





CONV 

90 

* 

INITIALIZATION AND BOUNDARY CONDITIONS 

I 



CONV 

91 

# 





CONV 

92 


IC . = ' 1 




CONV 

93 


DTDY ( 1 ) =0 »0 




. CONV 

94 


Y < 1 ) =0.0 




CONV 

95 


TRACT = 0.0 




CONV 

96 


FLASHsO. 0 




CONV 

97 


P I =3. 1 4 1 5927 i- ' 




CONV 

98 


DO 10 1=1 iNP 

, 



CONV 

99 


NEWT ( 1 1=1 • *< l.-YII ) >+TW 




CONV 

1 00 


IF CI.LT.NP) Y( 1+1 )=Y( I )+l./FLOAT(NH) 




•CONV 

10 1 


10 CONTINUE 




CONV 

102 

* 





CONV 

103 


R=R1*R2/(R1+R2) 




CONV 

1 04 


E = 2./( < 1 .-POI S1*P0IS1 )/El + ( 1 . — PO I S2*PO I S2 ) /E2 > 



CONV 

105 


B = 4 . *R/E*PHZ 




CONV 

106 


FL ASHK=0 . 24/SQRT ( P I *PHZ*DK*DRHO*DC*U*R/F ) 




CONV 

107 


TRANS! 1 ) =B/U/6. 




CONV 

108 


TRANS(Z) = B/U/6 . 




. CONV 

109 


TRANS! 3) =B/U/2. 




CONV 

1 10 


TRANS ! 4 ) = TRANS ( 3 ) 




CONV 

1 1 1 


TRANS ( B ) = TRANS ( 2 ) 




CONV 

1 12 


TRANS ! 6 ) sTRANS! 1 ) 




CONV 

1 13 


DELXB=B/FLOAT(NIP> 




CONV 

1 14 

* 





CONV 

1 15 

* 

CALCULATION OF LOAD 




CONV 

1 16 

-* 





CONV 

1 17 
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W=PHZ*PHZ*P I * ( 2 . /E ) *R 

CONV 

118 



PRINT 2 « PHZ «WtR*E 

conv 

1 19 


2 

FORMAT ( * 1 CONTROL PHZ = *E15.8»* W = *E15.8«* R = *E15.8, 

CONV 

120 



1* E = *E15.8) 

CONV 

121 

* 



CONV 

122 

* 


CALCULATION OF HALF-FILM THICKNESS 

CONV 

123 

* 



CONV 

124 



ET AENT = 1 0 • **< HERSA+HERSB*ALOG 1 0 ( TO I L“460 ♦ ) ) * 1 • 45E-7 

CONV 

125 



H- 1 .6*ALPHA**0.6*<ETAENT*U)**0«7*E**0,03*R**0»43/W*-*0. 13 

CONV 

126 



H= 1 . 2*H 

CONV 

127 



H=PHIT*H 

CONV 

128 



H=0.5*H 

CONV 

129 



print i « h» alpha. etaent.u.ch. ah 

CONV 

130 


1 

FORMAT ( #0H = *E12.5.* ALPHA = *E12»5»* ETAENT = *E12.5* 

CONV 

131 



1* U = *F6.0**CH = *FS.3«* AH = *F7»3> 

CONV 

132 

* 



CONV 

133 

* 


SLIP LOOP 

CONV 

134 

* 



CONV 

135 



NSL I P— 1 5 

CONV 

136 



DATA (SLIP! IU) . IU=1 *20)/»5. 1 . .2. .3. .4. »5» »6. #8. * 10». 15. *20. «30. 

•40CONV 

137 



1 . .50. *60 * ,5*0.0/ 

CONV 

138 



DO 6000 I U = 1 * NSL I P 

CONV 

139 



DATA ( TR ACTCF ( IP) . IP=1»20) / 2.0*0 ,0/ 

CONV 

140 



FLASH=FLASHK*TRACT*SLIP< IU) 

CONV 

141 



IF (FLASH.LT. 0.5) FLASH=0.0 

CONV 

142 



TW= TOI L+FLASH 

CONV 

143 



PRINT 6.FLASH ' 

CONV 

144 


6 

FORMAT ( * CONTROL FLASH = -»E15.8) 

CONV 

145 



U2U 1 =0 ,5*SL I P ( IU) 

CONV 

146 



PRINT 8. IU.SLIP< IU) 

CONV 

147 


8 

FORMAT (*1 CONTROL IU = *I3«* SLIP = -»E15.8) 

CONV 

148 

#■ 



CONV 

149 

* 


HERTZIAN PRESSURE LOOP 

CONV 

150 

•* 


: • • • - 

CONV 

151 



TRACT =0.0 

CONV 

152 



DO 5999 IP=1 .NIP2 

CONV 

153 



IPTRANS= IP 

CONV 

154 



XB= ( 2. -AFLOAT ( IP)-1 . > /£. /FLOAT INI P) 

CONV 

155 



P=PHZ*SQRT(XB*(2*-XB) ) 

CONV 

156 



DO 3 I = 1 . NP 

CONV 

157 



IF (IP.EO.l) T JM 1 ( I ) =TO I L 

CONV 

158 



IF (IP.GT.l) TJM1 ( I ) =NEWT< I) 

CONV 

159 


3 

CONTINUE 

CONV 

160 

* 



CONV 

161 

* 


SOLVE MOMENTUM EQUATION 

CONV 

162 

* 



CONV 

163 



I TCOUNT ~ 0 

CONV 

164 



IT = 0 

CONV 

165 



I TC= 1 

CONV 

166 



TM Ax=950 . = 

CONV 

167 



TCI) =95 1 . 

CONV 

168 


4 

IF (PRNT.EQ.2HON) PRINT 44 , I T 

CONV 

169 


44 

FORMAT (*0 IT = *13) 

CONV 

170 



IF (ITC.EO.O) GO TO 4970 

CONV 

171 



• IT= IT+1 

CONV 

172 



'IF (IT.GT.l) TMAX=T<IT-1) 

CONV 

173 



TMIn=TOIL 

CONV 

174 



I TC = 0 

CONV 

175 

4970 

DO 1 1 I = I T * NP 

CONV 

176 
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IF ( ITCOUNT.EQ.O) TEMT=NEWT(I) 


CONV 

177 


IF ( ITCOUNT.GT.O) TEMT=0 .5* ( T < I ) +NEWT ( I ) ) 


CONV 

178 


IF ( ITCOUNT.GT.200) GO TO 6003 


CONV 

179 


IF ( IP.LE.NIP.AND.TEMT.LT.TJM1 ( I ) ) 4980*4979 


CONV 

180 

4979 

IF (IP.GT.NIP. AND • TEMT « LT « T JM 1(11) 4982*4985 


CONV 

181 

4980 

TEMT=0.5*(TJM1 ( I )+T< I } ) 


CONV 

182 


GO TO 4985 


CONV 

183 

4982 

TEMT = 0.5* (TOIL+Yt I )*FLASH+Tt I 1 1 


CONV 

184 

4985 

IF (I.GT.IT) GO TO 4987 


CONV 

185 

* 

UPDATE TMAX AND TMlN 


CONV 

186 


IF ( TEMT.LT.TI IT} .AND. T{ IT} . LT . TMAX . AND. T ( I T ) .GT.TMIN) 

TMAX=T( IT) 

CONV 

187 


IF (TEMT.GT.TI IT} . AND . T { IT} . GT « TM I N . AND . T ( IT) . LT • TMAX ) 

TM I N=T ( IT) 

CONV 

188 


IF (PRNT.EQ.2HON) PRINT 4988 . TMAX * TM I N « TEMT 


CONV 

189 

4988 

FORMAT ( *OTMAX = *F9.4,* TM I N = *F9.4** TEMT = *F9. 

4) 

CONV 

190 

* 

BOUND TEMPERATURE 


CONV 

191 


IF (TEMT. GT. TMAX. OR • TEMT .LT . TM 1 N ) TEMT = 0 . 5* ( TMAX+TM I N ) 


CONV 

192 


IF (PRNT.EO.2HON) PRINT 4989 * TEMT . T ( I T > 


CONV 

193 

4989 

FORMAT!* TEMT = *F9.4«* T(IT) = *F9.4) 


CONV 

194 

* 

TEST FOR CONVERGENCE 


CONV 

195 


IF ( ABS( TEMT-T( IT) ) .LT. *3) 4986*4987 


CONV 

196 

4986 

I TC= 1 


CONV 

197 


IF (IT.EQ.NP) 5032.11 


CONV 

198 

4987 

T ( I ) =TEMT 


CONV 

199 


GINF( I > = 1 .2*P/( 2.52+. Ol 333* ( T ( I 1-492. ) >-l .45E4 


CONV 

200 


IF (GINF( I ) .LT.l .5 G INF ( I ) = 1 . 


CONV 

201 

* 

CALCULATION OF STEADY-STATE VISCOSITY 


CONV 

202 


ETEXP=ALPHA*P+(BETA+GAMMA*P)*(680.-T( I > }/'680.XT( I ) 


CONV 

203 


ETA2 ( I) = .62*EXP(ETEXP) *1 .45E-5 


CONV 

204 

* 

CALCULATION OF TIME-DEPENDENT VISCOSITY 


CONV 

205 


E T A 0 ( I ) = V I SC ( P. ETA2 ( I } . IVCODE ) 


CONV 

206 

1 1 

CONTINUE 


CONV 

207 


DUMMY=PSI 0(0.0) 


CONV 

208 


CALL SECANT(XK0*XK1 «TAU*PSI *XK*.00l *500* I TERR ) 


CONV 

209 


IF (ITERR.EQ.l) STOP 


CONV 

210 


IF (PRNT.EO.2HON) PRINT 41«TAU 


CONV 

21 1 

41 

FORMAT ( *OTAU = *E15.8) 


CONV 

212 

* 



CONV 

213 

* 

SOLVE ENERGY EQUATION 


CONV 

214 

* 



CONV 

215 


DO 5000 I = 1 »NP 


CONV 

216 


Q ( I ) =-TAU*DUDY ( I ) /COND*H 


CONV 

217 


QC=+OILRHO*OILC/COND*H*H*U*(T ( I l-TJMl ( I ) ) /DELXB 


CONV 

218 


Q< I ) =Q< I ) +QC 


CONV 

219 

5000 

CONTINUE 


CONV 

220 


IF (PRNT.EO.2HON) CALL PR1 


CONV 

221 


DO 5010 i=a,NP 


CONV 

222 


CALL INTEG1 ( 0. . Y( I ) « Y , Q * NP * DTD Y ( I > *C.GC,LLLL. I ERR) 


CONV 

223 

5010 

IF( IERR.NE.O) PRINT 5012* IERR 


CONV 

224 

5012 

FORMAT!* INTEGRATING Q. IERR = *13) 


CONV 

225 


DO 5020 1=2 *.NP 


CONV 

226 


CALL INTEG1 (0. , Y( I ) * Y , DTDY « NP * NEWT ( I ) « C * GC * LLLL * IERR) 


CONV 

227 

5020 

IF ( IERR.NE.O) PRINT 5022* IERR 


CONV 

228 

5022 

FORMAT ( * INTEGRATING DTDY. IERR = *13) 


CONV 

229 


XK7=TW— NEWT(NP) 


CONV 

230 


I TCOUNT= I TCOUNT+1 


CONV 

231 


NEWT ( 1 > = 0.0 


CONV 

232 


DO 5030 1=1, NP 


CONV 

233 

5030 

NEWT ( I ) = NEWT ( I )+XK7 


CONV 

234 
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GO TO 4 

5032 DO 5033 1=1 «NP 

5033 NEWTC I ) =T ( 1) 

IF (IVCODE.NE.O) PRINT 5034 

5034 F0RmAT(*0 MINIMUM VISCOSITY REDUCTION*) 

CALL PR 1 

PRINT 99 

PRINT 5040. < (NEWT! I ) ♦ I ) • 1 = 1 *NP) 

5040 FORMAT < * NEW TEMP = *F7.2,* I = *14) 

* 

* CALCULATION OF TRACTION 

* 

TR ACT = TRACT+8 /FLOAT (NIP) *T AU 

5999 CONTINUE 

TRACTCF ( I U) =TRACT/W 
PRINT 600 2 *TRACTCF ( IU) 

6002 FORMAT (*OCONTROL TRACTCF(IU) = *E15»8) 

# 

* PLOT TEMPERATURE PROFILE 

# 

IF (nGRAPH.EQ.O) GO TO 6000 
DO 5052 NGR=1.NGRAPH 
PRINT 50 5 G * SL I P ( I U ) 

5050 FORMAT(*l(U2 - U1 ) = *E15.8) 

PRINT 99 

CALL STPLT1 ( 1 .NEWT.Y.NP. 1H** 1 • 1 HY > 

PRINT 5051 

5051 FORMAT ( 1 HC .92X, 1 1 HTEMPER ATURE > 

5052 CONTINUE 

6000 CONTINUE 

6003 PRINT 98 

PRINT 60 0 1 * < ( TRACTCF (IU)*SLIP(IU))*IU=1 *NSLIP) 

6001 FORMAT ( * CONTROL TRACT = *E15.8** U2-U1 = *E15»8> 

PUNCH 6005* I (SLIP! IU ) * TRACTCF ( IU) ) • IU=1 «NSL1P) 

6005 FORMAT ( 2F 10*5) 

■*, 

* PLOT TRACTION COEFFICIENT 

* 


IF (mGRAPH.EO.O) GO TO 6013 
DO 6012 MGR=1*MGRAPH 
PRINT 98 

CALL STPLT1 ( 1 .SLIP. TRACTCF, NSLIP* 1H** 14, 14HTRACTI0N COEF.) 
PRINT 6011 

60 1 1 FORMAT < 1 HO * 1 OOX , 7HU2 - Ul) 

6012 CONTINUE 

6013 CONTINUE 

98 FORMAT (*1*) 

99 FORMAT <*0*) 

GO TO 9000 

9999 CONTINUE 
STOP 
END . 
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CONV 245 
CONV 246 
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CONV 256 
CONV 257 
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CONV 277 
CONV 278 
CONV 279 
CONV 280 
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FUNCTION PS I (TAU) 



PS IBL 

1 

* 

-* * 

-fr-***--*-****-**-*-**-***-**# 

* 

* * * 

*-***-***** PS |BL 

2 

* 





*PS IBL 

3 

* 


FUNCTION PS I 



*PS I BL 

4 

-H- 





*PS I BL 

5 

* 





*PS I BL 

6 

-H- 


THIS FUNCTION SUBROUTINE SUPPLIES THE RHEOLOGICAL MODEL *PSIBL 

7 

■K- 


FOR THE LUBRICANT. IN THIS CASE* WE ARE 

USING THE *PS I BL 

8 

* 


BARLOW AND LAMB VISCOELASTIC MODEL WITH 

A 

LIMITING SHEAR STRESS. *P S I BL 

9 

* 





*PS IBL 

10 

* 

* * 

*-***-«•*■*-**-*#**#*#■*•*■** 

* 

* * * 

*****•*##•* PS IBL 

1 1 



COMMON C.GC.LLLL 



PS IBL 

12 



COMMON /CPSI/ GINF*ETA0«T.Y*H«U2U1 *NP«DUDV 

, IC*TG 

, OMEGA, CH* AH PS IBL 

13 



COMMON /CPS I 0/ XK l ♦ XKO 



PS IBL 

14 



DIMENSION G I NF (21) *ETA&<21 )*T(21).Y(21>, 

DUDY (21 ) 

*TG(21 ) .OMEGA (21 ) PS IBL 

15 



DIMENSION C(2*3*21 ) » GC (21*21 ) 



PS I BL 

16 



real intg 



PS IBL 

17 



PRNT=3HOFF 



PS IBL 

18 


100 

CONTINUE 



PS IBL 

19 



IF (PRNT.EQ.2H0N) PRINT 200.TAU 



PS IBL 

20 


200 

FORMAT(*OTAU = *E15.8) 



PS IBL 

21 



DO 210 1 = 1* NP 



PS IBL 

22 



TG ( I >=TAU/GINF( I ) 



PS IBL 

23 



IF <TG( I ) .LT.0.25) GO TO 205 



PS IBL 

24 

* 


SLIP REGION 



PS IBL 

25 



OMEGA ( I ) = ( TG ( I ) -.24999631*1 .E6 



PS IBL 

26 



GO TO 208 



: PS IBL 

27 

* 


VISCOELASTIC REGION ' 



PS IBL 

28 


205 

OMEGA ( I )=55.2*TG( I )*TG< I >+TG( I ) 



PS IBL 

29 


208 

DUDY ( I ) = G I NF ( I ) /E T A 0 ( I >*OMEGA( I )*H 



PS IBL 

30 


210 

CONTINUE 



PS IBL 

31 



GO TO (220.230) . IC 



PS I BL 

32 


220 

CALL INTEG ( 0 . , 1 . * Y ♦ DUD Y * NP , I NTG , C , GC * LLLL 

, IERR) 

PS IBL 

33 



IC = 2 



PS I BL 

34 



GO TO 240 



PS I BL 

35 


230 

CALL I NTEG2 ( 0 • . 1 . *Y*DUDY*NP, I NTG * C * GC * LLLL 

, IERR) 

PS IBL 

36 


240 

IF (IERR.NE.O) PRINT 24 1 « I ERR . TAU ' 



PS IBL 

37 


241 

FORMAT!* IERR = *14** AT TAU = *E15.8> 



’ PS I BL 

38 



PS I = ALOG 10(1 NTG/U2U1 ) 



PS IBL 

39 



IF (PRNT.EQ.2HON) PRINT 250* I NTG* PS I 



PS IBL 

40 


250 

FORMAT!* INTG = *E15.8.* PS I = *E15.8) 



PS IBL 

4 1 



TDUD Y = 0 . 0 



PS IBL 

42 



DO 310 1 = 1* NP 



PS I BL 

43 


310 

TDUDY=TDUDY+DUDY( I ) 



PS IBL 

44 



AD J = U2U 1 * FLO AT ( NP ) /TDUDY 



PS IBL 

45 



IF ( ADJ.GT.0.99) GO TO 390 



PS IBL 

46 



DO 320 1=1, NP 



PS IBL 

47 


320 

DUDY ( I > = AD J*DUDY ( I ) 



PS IBL 

48 


390 

CONTINUE 



PS IBL 

49 



RETURN 



PS IBL 

50 



ENTRY PS I 0 



PS IBL 

51 

* 


CHECK OMEGA = 3.7 



PS IBL 

52 



XK 1 =T AU= 0 . 25*G I NF ( 1 > 



PS IBL 

53 



DO 4 10 I = 1 *NP 



' - PS IBL 

54 



TG ( I ) =TAU/GINF( I ) 



PS IBL 

55 



OMEGA ( I ) =55 • 2*T G ( I ) *TG ( I )+TG( I ) 



PS IBL 

56 


4 1 0 

DUDY ( I ) =G INF ( I ) /ETAO ( I ) *OMEGA ( I ) *H 



PS IBL 

57 



CALL INTEG ( 0 . , 1 . * Y * DUD Y * NP « I NTG * C * GC * LLLL 

, IERR) 

PS IBL 

58 
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IF ( I ERR • NE * 0 > PRINT 241.IERR.TAU 

PS I8L 

59 

IF ( ALOG 10(1 NTG/U2U1 > ) 20.20.30 

PSIBL 

60 

SLIP MODEL 

PS IBL 

61 

XKO= 10.-*-XKl 

PSIBL 

62 

PSI =0.0 

PSIBL 

63 

RETURN 

. PSIBL 

64 

VISCOELASTIC MODEL 

PSIBL 

65 

XKO=0. l*XKl 

PSIBL 

66 

PS I =0.0 

PSIBL 

67 

RETURN 

PSIBL 

68 

END 

PSIBL 

69 





FUNCTION PSI(TAU) 



PS I MX 

1 

* 

* * 

********************** 

* 

* * 

******** PS I MX 

2 

* 





*PS I MX 

3 

* 


FUNCTION PS I 



*PS I MX 

4 

* 





♦PS I MX 

5 

* 





♦PS I MX 

6 

* 


THIS FUNCTION SUBROUTINE SUPPLIES THE RHEOLOGICAL 

MODEL *PS I MX 

7 

* 


FOR THE LUBRICANT* IN THIS CASE* WE ARE USING 

THE *PSIMX 

8 

* 


MAXWELL VISCOELASTIC MODEL WITH A LIMITING 

SHEAR 

STRESS* ♦PS I MX 

9 

* 





♦PS I MX 

10 

* 

* * 

********************** 

* 

* * 

♦ *♦*♦*** PS 1 MX 

1 1 



COMMON C«GC«LLLL 



PS IMX 

12 



COMMON /CPSI/ GINF«ETA0,T*Y«H*U2U1 *NP*DUDY* IC* 

TG, 

OMEGA, CH* AH PS I MX 

13 



COMMON /CPS 10/ XK 1 * XKO 



PS IMX 

14 



DIMENSION G I NF (21) *ETA0(21 )*T(21),YC21) «DUDY(21 > , 

TG< 2 1 > * OMEGA (21) PSIMX 

15 



DIMENSION C<2«3*21)*GC <21*21) 



PS I MX 

16 



REAL intg 



PSIMX 

17 



PRNT=3HOFF 



PSIMX 

18 


100 

CONTINUE 



PSIMX 

19 



IF ( PRNT • EQ* 2HON ) PRINT 200*TAU 



PSIMX 

20 


200 

FORMAT<*OTAU = *E15*8) 



PSIMX 

21 



DO 210 1 = 1* NP 



PSIMX 

22 



TG ( I ) =T AU/G I NF ( I ) 



PSIMX 

23 



IF (TG< I ) .LT.0.50) GO TO 205 



PSIMX 

24 

* 


SLIP REGION 



PSIMX 

25 



OMEGA ( I ) = ( TG t I ) -0*4999990 )♦ 1 *E6 



PSIMX 

26 



GO TO 208 



PSIMX 

27 

* 


VISCOELASTIC REGION 



PSIMX 

28 


205 

G I NFK = 0 • 5/TG < I ) 



PSIMX 

29 



OMEGA ( I 1 = 1./<GINFK+SQRT<GINFK*GINFK-I • ) ) 



PSIMX 

30 


208 

DUD Y ( I ) =G I NF < I )/ETAO< 1 )*OMEGA< I )*H 



PSIMX 

31 


210 

CONTINUE 



PSIMX 

32 



GO TO (220,230) , IC 



PSIMX 

33 


220 

CALL INTEG < 0 . , 1 * . Y * DUDY *NP * I NTG * C ♦ GC * LLLL * 

I ERR) 

PSIMX 

34 



(\J 

11 

u 
*— « 



PSIMX 

35 



GO TO 240 



PSIMX 

36 


230 

CALL I NTEG2 ( 0 • , 1 . *Y«DuDy*NP, I NTG ♦ C * GC , LLLL , 

I ERR) 

PSIMX 

37 


240 

IF ( I ERR • NE • 0 ) PRINT 241*IERR«TAU 



PSIMX 

38 


241 

FORMAT < * I ERR = *14** AT TAU = *E15*8) 



PSIMX 

39 



PS I = ALOG 1 0 ( I NTG/U2U 1 ) 



PSIMX 

40 



IF ( PRNT • EQ • 2HON ) PRINT 250*INTG*PSI 



PSIMX 

41 


250 

FORMATL* INTG = *E15.8,* PS I = *E15.8) 



PSIMX 

42 



TDUDY=0* 0 



PSIMX 

43 



DO 310 1=1, NP 



PSIMX 

44 


310 

TOUDY=TDUOY+DUDY< I ) 



PSIMX 

45 



AO J=U2U1 *FLOAT < NP ) /TDUDY 



PSIMX 

46 



IF < ADJ.GT.0.99) GO TO 390 



PS IMX 

47 



DO 320 1=1, NP 



PSIMX 

48 


320 

DUDY ( I ) =ADJ*DUDY( I ) 



PSIMX 

49 


390 

CONTINUE 



PSIMX 

50 



RETURN 



PSIMX 

51 



ENTRY PS 10 



PSIMX 

52 

* 


CHECK OMEGA =1 



PSIMX 

53 



XK 1 =TAU= 0 • 50*G I NF ( 1 ) 



PSIMX 

54 



DO 410 1=1, NP 



PSIMX 

55 



TG ( I )=TAU/GINF< I ) 



PSIMX 

56 



G I NFK = 0 » 5/TG ( I ) 



PSIMX 

57 



OMEGA < I ) = 1 ./(GINFK+SQRT(GINFK*GINFK-1 . ) ) 



PSIMX 

58 


198 



410 DUDY ( I >=GlNF< I )/ETA0< I >*OMEGA< I )*H 

CALL INTEG (O. i 1 . iYi DuDY»NP • 1 NTG »C ♦ GC *LLLL • 1 ERR ) 
IF ( IERR «NE« 0 ) PRINT 241«IERRtTAU 
IF ( ALOG 10(1 NTG/U2U1 ) ) 20«20«3 0 

* SLIP MODEL 
20 XK0=10»*XK1 

PSI=0»0 

RETURN 

* VISCOELASTIC MODEL 
30 XK0=0«1*XK1 

PSIaO.O 

RETURN 

END 


PS I MX 59 
PS I MX 60 
PS 1 MX 61 
PS1MX 62 
PS I MX 63 
PS I MX 64 
PS I MX 65 
PS I MX 66 
PS I MX 67 
PS I MX 68 
PS I MX 69 
PS I MX 70 
PS I MX 71 


199 



SUBROUTINE RTMI<X, F. FCT • XL I » XR I * EPS* IEND, IER) 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 
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c 

c 

c 
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c 
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SUBROUTINE RTMI 
PURPOSE 

TO SOLVE GENERAL NONLINEAR EQUATIONS OF THE FORM FCT(X)=0 
BY MEANS OF MULLER-S ITERATION METHOD* 

USAGE 

CALL RTMI <X*F*FCT»XLI » XR I lEPS. IEND* IER) 

PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT* 

DESCRIPTION OF PARAMETERS 

X - RESULTANT ROOT OF EQUATION FCT(X>=0. 

F - RESULTANT FUNCTION VALUE AT ROOT X. 

FCT - NAME OF THE EXTERNAL FUNCTION SUBPROGRAM USED. 

XL I - INPUT VALUE WHICH SPECIFIES THE INITIAL LEFT BOUND 

OF THE ROOT X. 

XR I - INPUT VALUE WHICH SPECIFIES THE INITIAL RIGHT BOUND 
OF THE ROOT X. 

EPS - INPUT VALUE WHICH SPECIFIES THE UPPER BOUND OF THE 
ERROR OF RESULT X* 

IEND - MAXIMUM NUMBER OF ITERATION STEPS SPECIFIED. 

IER - RESULTANT ERROR PARAMETER CODED AS FOLLOWS 
IER=0 - NO ERROR* 

I ER= 1 - NO CONVERGENCE AFTER IEND ITERATION STEPS 
FOLLOWED BY IEND SUCCESSIVE STEPS OF 
BISECTION. 

IER=2 - BASIC ASSUMPTION FCT ( XL I ) *FCT ( XR I > LESS 
THAN OR EQUAL TO ZERO IS NOT SATISFIED* 

REMARKS 

THE PROCEDURE ASSUMES THAT FUNCTION VALUES AT INITIAL 
BOUNDS XL I AND XR I HAVE NOT THE SAME SIGN. IF THIS BASIC 
ASSUMPTION IS NOT SATISFIED BY INPUT VALUES XL I AND XR I • THE 
PROCEDURE IS BYPASSED AND GIVES THE ERROR MESSAGE IER=2. 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 

THE EXTERNAL FUNCTION SUBPROGRAM FCT(X) MUST BE FURNISHED 
BY THE USER. 

METHOD 

SOLUTION OF EQUATION FCT(X)=0 IS DONE BY MEANS OF MULLER-S 
ITERATION METHOD OF SUCCESSIVE BISECTIONS AND INVERSE 
PARABOLIC INTERPOLATION, WHICH STAPTS AT THE INITIAL BOUNDS 
XL I AND XRI. CONVERGENCE IS QUADRATIC IF THE DERIVATIVE OF 
FCT < X ) AT ROOT X IS NOT EQUAL TO ZERO. ONE ITERATION STEP 
REQUIRES TWO EVALUATIONS OF FCT(X), FOR TEST ON SATISFACTORY 
ACCURACY SEE FORMULAE (3*4) OF MATHEMATICAL DESCRIPTION. 

FOR REFERENCE, SEE G. K. KRISTIANSEN* ZERO OF ARBITRARY 
FUNCTION, BIT, VOL . 3 (1963), PP. 205-206. 
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c 

c 

c 

C PREPARE ITERATION 

IER = 0 
XL=XLI 
XR=XRI 
X = XL 
TOL = X 
F = ECT< TOL ) 

IF<F) 1 « 16« 1 

1 fl=f 

X = XR 
TOL=X 
F=FCT(TOLl 
IF(F>2,16.2 

2 FR=F 

IF (SIGN! 1 « .FLl+S IGN( 1..FRI ) 25 ,3,25 

BASIC ASSUMPTION FL*FR LESS THAN 0 IS SATISFIED. 
GENERATE TOLERANCE FOR FUNCTION VALUES. 

3 I =0 

TOLF = 1 00 . *EPS 


START ITERATION LOOP 
4 1 = 1+1 

START BISECTION LOOP 
DO 13 K=U I END 
X=.5*tXL+XR> 

TOL = X 
F = FCT ( TOL 1 
IF ( F ) 5 » 1 6 « 5 

8 IF (SIGN! 1 • ,F)+S IGN< 1 . ,FR) >7.6.7 

INTERCHANGE XL AND XR IN ORDER TO GET THE SAME SIGN IN F AND FR 

6 TOL =XL 
XL = XR 
XR=TOL 
TOL=FL 
FL='FR 
FR = TOL 

7 TOL=F-FL 
A=F*TOL 
A = A + A 

IF (A-FRtt(FR-FL) 18.9.9 

8 IF ( I-IEND) 17, 17,9 

9 XR = X 
FR = F 

TEST ON SATISFACTORY ACCURACY IN BISECTION LOOP 
TOL = EPS 
A= ABS ( XR ) 

IF ( A-l • ) 11,11,10 

10 TOL=TOL*A 

11 IF(ABS(XR-XL)-TOL) 12, 12. 13 

12 IFC A.BS(FR-FL)-TOLF) 1^, 14, 13 
•13 CONTINUE 
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C END OF BISECTION LOOP 

C 

C NO convergence after i end iteration steps followed by I end 

C SUCCESSIVE STEPS OF BISECTION OR STEADILY INCREASING FUNCTION 

C VALUES AT RIGHT BOUNDS. ERROR RETURN. 

I ER = 1 

14 IF(ABS(FR)-ABS(FL) ) 16, 16. 15 

15 X=XL 

f=fl 

16 RETURN 

COMPUTATION OF ITERATED X-VALUE BY INVERSE PARABOLIC INTERPOLATION 

17 A = F R— F 

DX= (X-XL ) *FL*< 1 . +F * ( A — TOL ) / ( A* ( FR~ FL ) ) > /TOL 
XM = X 
FM = F 
x=xl-dx 
T OL = X 

F=FCT( TOL ) 

IF(F) 18, 16. 18 

TEST ON SATISFACTORY ACCURACY IN ITERATION LOOP 

18 TOL=EPS 
A = ABS < X ) 

IF ( A-l • ) 20,20. 19 

19 TOL=TOL*A 

20 IF( ABS(DX)-T0L)21 ,21 ,22 

21 IF(ABS(F)-TOLF) 16. 16,22 

PREPARATION OF NEXT BISECTION LOOP 

22 IF ( SIGNC 1 ♦ ,F )+SIGN( 1 . , FL > ) 24 . 23 , 24 

23 XR=X 
FR = F 
GO TO 4 

24 XL=X 

fl=f 

XR=XM 
FR = FM 
GO TO 4 
C END OF ITERATION LOOP 

C 
C 

C ERROR RETURN IN CASE OF WRONG INPUT DATA 

25 IER=2 
RETURN 
END 
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APPENDIX C 


NOMENCLATURE 

a rise parameter of the hyperbolic model 

b = 4P Hz R/ E > half -width of Hertzian contact 

c limiting shear stress/limiting shear modulus ratio 

c specific heat of the lubricant 

c^ specific heat of the disk 

D shear rate 

2 2 

E = 5f(l-v^/ E ^) + effective modulus of elasticity 

of the disks 

E^jE^ elastic moduli of the two disks 

Ei(x) exponential integral 

f fractional free volume 

£2 equilibrium free volume 

* 

G complex shear modulus 

G high frequency limiting shear modulus 

00 

G = G /K 

CO CO 

h ha If -para llel lubricant film thickness 

h minimum lubricant film thickness 

o 

k thermal conductivity of the lubricant 

k thermal conductivity of the disk 

m 

K Oldroyd-Dyson parameter 

K bulk modulus 

low frequency bulk modulus 

K high frequency bulk modulus 

00 

bulk modulus associated with molecular rearrangements in 
free volume 

complex relaxational modulus 

203 


high frequency value of 

torsional spring constant 

J inertia 

* 

J complex compliance 

p pressure 

p^ z maximum Hertzian pressure 

P maximum Hertzian pressure (on graphs) 

P pressure step 

P normal stress 

R = RjR^/CR^ + Rj), effective radius of the disk pair 

R^,R 2 radii of the disks 

s = ln(T] 2 /7]) 

s x = ln(Tl 2 /Tl 1 ) 

t time 

T temperature 

T lubricant inlet temperature (on graphs) 

T, bulk temperature of the disk 

b 

T ent lubricant entrance temperature 

T s mean surface temperature of the disk in the contact zone, 
"flash temperature" 

T q reference temperature at which there is no free volume 
u,v velocity components in the fluid film 
U mean rolling speed 

U^,U 2 surface speeds of the disks 

'k — 

U reference rolling speed at which T] = \ T\j_q 

v specific volume 

v^ free volume 
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V. 

1 


w 


x,y,z 

* 

Z 

a 


Y 

Y 
T1 

\nt 

If 

\ 

% 

Tl 2 

ii 

x 


p- 

* 

p* 

V l» v 2 

§ 

P 


specific volume after elastic deformation only 

occupied volume 

initial specific volume 

equilibrium specific volume 

load per unit length of cylinder 

Cartesian coordinates 

complex shear mechanical impedance 

viscosity-pressure coefficient 

shear strain 

shear rate 

shear viscosity 

viscosity of the lubricant at entrance conditions 

free volume viscosity 

volume viscosity 

initial viscosity 

equilibrium viscosity 

effective viscosity 

Maxwell relaxation time 

retardation time 

volume relaxation time 

coefficient of friction 

complex fluidity 

Poisson's ratio for the two disks 
slide/roll ratio 
density of the lubricant 
density of the disk 
shear stress 
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film thickness thermal reduction factor 
angular frequency 


non-dimensional shear rate 
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